{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "aca796a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import pymc as pm\n",
    "import arviz as az\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "from statsmodels.stats.multitest import multipletests\n",
    "from scipy.stats import rankdata\n",
    "import matplotlib.patches as mpatches\n",
    "import glob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "7179ca98",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the comparison you want to analyse\n",
    "# Example: Comparison (b): cdd_no_fft vs basic_no_fft\n",
    "model_1 = \"basic_no_fft\"\n",
    "model_2 = \"basic_fft\"\n",
    "metric = \"Specificity\"  # Change to any other metric (Accuracy, Precision, etc.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "9a86a2c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load your CSV\n",
    "df = pd.read_csv(\"results_metrics.csv\")\n",
    "\n",
    "# Filter to External dataset\n",
    "df = df[df[\"Dataset\"] == \"External\"]\n",
    "\n",
    "# Standardise model names\n",
    "df[\"Model\"] = df[\"Model\"].replace({\n",
    "    \"Basic Model No FFT\": \"basic_no_fft\",\n",
    "    \"Basic Model With FFT\": \"basic_fft\",\n",
    "    \"CDD-based Model No FFT\": \"cdd_no_fft\",\n",
    "    \"CDD-based Model With FFT\": \"cdd_fft\"\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0c9ec035",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2e2776aa5bef4c1eb6b19816048215ce",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 23 seconds.\n"
     ]
    }
   ],
   "source": [
    "#GLM-style Bayesian model for two groups\n",
    "\n",
    "# Extract the metric values for the two models\n",
    "y1 = df[df[\"Model\"] == model_1][metric].values\n",
    "y2 = df[df[\"Model\"] == model_2][metric].values\n",
    "\n",
    "\n",
    "# Combine data for both groups\n",
    "y = np.concatenate([y1, y2])\n",
    "group = np.concatenate([\n",
    "    np.zeros_like(y1),  # 0 = baseline (e.g. basic_no_fft)\n",
    "    np.ones_like(y2)    # 1 = treatment (e.g. cdd_no_fft)\n",
    "])\n",
    "\n",
    "# Build GLM model\n",
    "with pm.Model() as model:\n",
    "    # Intercept: baseline mean\n",
    "    beta0 = pm.Normal(\"beta0\", mu=0.5, sigma=1)\n",
    "\n",
    "    # Group effect (difference between group 1 and group 0)\n",
    "    beta1 = pm.Normal(\"beta1\", mu=0.0, sigma=1)\n",
    "\n",
    "    # Shared standard deviation\n",
    "    sigma = pm.Exponential(\"sigma\", lam=10)\n",
    "\n",
    "    # Linear predictor\n",
    "    mu = beta0 + beta1 * group\n",
    "\n",
    "    # Likelihood\n",
    "    y_obs = pm.Normal(\"y_obs\", mu=mu, sigma=sigma, observed=y)\n",
    "\n",
    "    # The treatment effect is directly beta1\n",
    "    trace = pm.sample(2000, tune=1000, chains=4, target_accept=0.95, return_inferencedata=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3b2a4907",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.microsoft.datawrangler.viewer.v0+json": {
       "columns": [
        {
         "name": "index",
         "rawType": "object",
         "type": "string"
        },
        {
         "name": "mean",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "sd",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "hdi_3%",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "hdi_97%",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "mcse_mean",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "mcse_sd",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "ess_bulk",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "ess_tail",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "r_hat",
         "rawType": "float64",
         "type": "float"
        }
       ],
       "ref": "ff2cbc70-9094-49d9-b30c-7e9b7b7c4daa",
       "rows": [
        [
         "beta0",
         "0.925",
         "0.003",
         "0.921",
         "0.931",
         "0.0",
         "0.0",
         "3252.0",
         "3727.0",
         "1.0"
        ],
        [
         "beta1",
         "0.033",
         "0.004",
         "0.026",
         "0.04",
         "0.0",
         "0.0",
         "3371.0",
         "4043.0",
         "1.0"
        ],
        [
         "sigma",
         "0.008",
         "0.002",
         "0.006",
         "0.011",
         "0.0",
         "0.0",
         "4811.0",
         "4758.0",
         "1.0"
        ]
       ],
       "shape": {
        "columns": 9,
        "rows": 3
       }
      },
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>beta0</th>\n",
       "      <td>0.925</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.921</td>\n",
       "      <td>0.931</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3252.0</td>\n",
       "      <td>3727.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta1</th>\n",
       "      <td>0.033</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.026</td>\n",
       "      <td>0.040</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3371.0</td>\n",
       "      <td>4043.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>0.008</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.011</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4811.0</td>\n",
       "      <td>4758.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\n",
       "beta0  0.925  0.003   0.921    0.931        0.0      0.0    3252.0    3727.0   \n",
       "beta1  0.033  0.004   0.026    0.040        0.0      0.0    3371.0    4043.0   \n",
       "sigma  0.008  0.002   0.006    0.011        0.0      0.0    4811.0    4758.0   \n",
       "\n",
       "       r_hat  \n",
       "beta0    1.0  \n",
       "beta1    1.0  \n",
       "sigma    1.0  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(trace, var_names=[\"beta0\", \"beta1\", \"sigma\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "00902dd6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAAG3CAYAAACjR8/oAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbNNJREFUeJzt3Qd4FNXaB/A3vRfS6IEk9A4ivQkCimIFRVGx3Iu9cy2f116uvXv1Xq9ibwh2BBRReu+9JUACIaGk92S+5382s2w2u6mbzJb/73nCLrO7s2dnZ2bfOec953hpmqYJERERkQvzNroARERERI3FgIaIiIhcHgMaIiIicnkMaIiIiMjlMaAhIiIil8eAhoiIiFweAxoiIiJyeQxoiIiIyOUxoCEiIiKX12QBTV5envztb3+TVq1aiZeXl9xzzz1q+fHjx2XKlCkSHR2tlr/++utNVQSqp44dO8r1119fZdm+fftkwoQJEhERob6v77//Xj766CN1PyUlpV7rx7rxHs3lxRdflG7duklFRYUYqT6f+4knnlDb1pXpn+HEiRPN+r54T7y3o5WVlckDDzwg7du3F29vb7nkkktqPMfVx5gxY6RXr17iaXDuwDZ7+eWXxVkY8V38+eefajvg1h0ct/P7but3xJ6dO3eKr6+vbN++vWkDGv2HzN7f6tWrzc997rnn1PNvvfVW+fTTT+Xaa69Vy++9915ZuHChPPzww2r5eeedJ46G965pgxl5EN9www2SlJQkgYGB6kQ4atQoefzxx8VZzZgxQ7Zt2ybPPvus+r4GDhzosHUXFBSoH6CmOJhzcnLkhRdekAcffFD9CDmTpvzc5HgffvihvPTSS+pE/fHHH6tzmL1z3MqVK9V3m5WVZVh57Z3/jCjb/PnzmyTIJOd0r53fd1u/I1988YXNCo0ePXrIBRdcII899lj9C6DVw+zZszHvk/bUU09pn376abW/zMxM83MHDx6sDR8+vNo6WrZsqU2fPl1rSiEhIdqMGTM0Z7Jv3z4tMjJSa926tfbII49o77//vtqOl1xyiRYQEKA5g6KiIq2kpMT8/4KCAvV9o7yWysrKtMLCQq2ioqJe68e68R467C9Y/+OPP6452muvvaaFh4erchqtPp+7tLTUKcrcGPhc+HyW54PmgO2G7edoV155pda2bdtqy22d41566SX12ZOTk+u07tGjR2s9e/bUmuP8V9+yOcLtt9+u3tMayoDlKJOzaIrvojbl5eVqv8WtO2hp4/fd3u/IBRdcoHXo0MHmeubPn69es3///nq9v29DorDzzz+/1iv1jIwMFWnZWh4ZGSme5rXXXlNV1Js3b5YOHTpU2ybOICAgoMr/MzMz1a319+Xj46P+6svPz0+ay+zZs+Wiiy5SNWFGq8/nRlUr/qj+muq7tnfOsneO8zSocQwODja6GC4JtcfOcI5qymPF3u9ITc4991xp0aKFqhF96qmn6l6AhtTQrFu3zu5zlixZop5j/ae/1vpPd/r0ae3uu+/W2rVrp/n7+2tJSUna888/Xy1yxf9ff/11rVevXqpmIyYmRps4caK5TLbeo7bamuPHj2s33nijFhcXp9bZp08f7aOPPrJ7RfGf//xHS0xMVOUcOHCgtnbt2lq3HcrYsWNHrS4QtSJ6Xbhwoda3b19Vpu7du2tz586t9lxHbTf9ffVtpV9lW/7p0bT+XVpf6SGqHjVqlBYaGqqFhYWpbfP555+bH8e69XXo29P6D+/74YcfqvsbN26s9nmfffZZzdvbW0tNTbW7/Q4ePKheX9N3+Oqrr2rx8fFaYGCgKvO2bduqrWfx4sXaiBEjtODgYC0iIkK76KKLtJ07d1Z5Tk5Ojtr++FzY/rGxsdq5556rbdiwod6f23K763DFOGbMmGplw/fZpk0b7fLLL6+yDDVTPXr0UN8x9ueZM2dqp06d0pqT/hl27dqlTZ06Ve0LUVFR2l133VWt9gnf9TnnnKO2G7Yf9vN///vf1daJ/XTChAladHS0+s5wLN1www1VnmOr1gv7CY5t1Ixi/XjdLbfcohUXF9f6Oex9V/bOcfiebS2vqUZErxVYv369NnToUPNne/fdd6s9F7V8jz32mDrG8VlwzP/jH/+oUvtnr1y2jmfrsqGWfcCAAaoMLVq0UDVThw8ftlvekSNHakFBQWr/t8Xe9qjv+XTLli1qXQkJCWq/Ri0AvvsTJ05UeZ7+GVEbjufjmEUt7fXXX6/l5+fX+n3X9bvAvvPoo4+qbYX14/yA88Qff/xRbZ1ffvmlep5+TsT5F+dhnb4v4dbS6tWrtfPPP1/V6mP9vXv3rvK62iypXO/XX3+tPfPMM6qGEdtu7NixavtY++abb8zfPY4x1LLUdI61Zu/33d7vCLa1vd8X3aWXXqp+i+ujQZeC2dnZ1RL+kEODRKDu3burNjK0pbVr107uv/9+9Xj//v3N7czjx4+X6667rkqEP3r0aElLS5Obb75Z4uPjVXsv2uGOHTtWpZ3tpptuUu3WqCVCQh4S9pYtW6byd1BrhPfA8kGDBsnMmTPVa5CzYk9hYaFKBtu/f7/ccccdkpCQIHPmzFGJnGhrvvvuu6s8H+1+ubm5qpz4zEg8veyyy+TgwYM1XomjVub333+XP/74Q8aOHVvrNkYS1ZVXXim33HKLan9EjcPUqVNlwYIFavs5ertZw2dCRI3v8aqrrpJJkyZJaGio3fJi3TfeeKP07NlTvT9eu2nTJlXeq6++utrzY2Nj5d1331X5B5deeql6P+jTp4/6Dm6//Xb5/PPP1X5jCcvwfbVt29ZuWbANYMCAATYf/+STT9R3iPcoKiqSN954Q30naONt2bKleg6+K2yrxMRElQOA/eStt96S4cOHy8aNG81Jvvh+vv32W7Xv4Gr95MmTsnz5ctm1a5fN96/pc9uCfQDvn56ernKudHiPo0ePyrRp08zLsA/ge0Ce1l133SXJycny9ttvq+9hxYoVNe6fxcXFapvURUxMTJ2ed8UVV6jt9K9//UvtZ2+++aacPn1abX8dtgX2GdSmoWbqp59+kttuu00lcuP70a/6kFCIbffQQw+pfQv5aPPmzavx/bF9cB7AcYxzARLEcazg+8Kx4+/vX+Pr8X44n6DdH7Wr+Bxg7xzXu3dvKSkpkS+//FLVyOrbCeupCbYJji9sLxxr33zzjdo/UD4cU4DtgW2E7x2fBWXA/or32bt3rzlnxt75LyQkRD3PXtnwGR999FFVBrweV9XY35Hjh/3H8uoa+ziODex711xzjfmYsYb9Ed/Bb7/9psplS13Op3g9/o/9GsfAjh075L///a+6xX5lnUSPz4BzCL4vHKv/+9//JC4uTuXU1aYu3wXy87BOPP73v/9dlf+DDz6QiRMnytq1a6Vfv37mcuM548aNM783zgs4Fq1/VyzhdRdeeKG0bt1aPQ+fGa/7+eefa3ydLc8//7yqBZo1a5b63cb2nT59uqxZs8b8HP2ccfbZZ6tthsRenBNRTuvv3h7sJ7Z+33Fes/U7gv0R5UlNTVX7I1j/vpx11lnyww8/qO0dHh4udVKf6MdeFIY/6zwQvZbBGp6LdlVLTz/9tGr33bt3b5XlDz30kObj42O+SkAEjNfjSs+aZT5HfXJoEPVinZ999lmVnAdE6IiqcQVueUWB6NXyiveHH35Qy3/66aca32f79u3qagbP7devn7qq+f77721eOWDb4XmWNTLZ2dnqKrN///5Ntt0sa2hqaue2rqHJyspSVx/IKbC+Ardcv2VNRW25JFdddZWqgbCsaUKNjV7bV5N//vOf6nm5ublVluufB9+D5dXHmjVr1PJ7773XvAzfEWo4Tp48WeVKEbVD1113nXkZrgKt92dr9fnc1jU0e/bsUf9/6623qjzvtttuU/sn2qdh2bJl6nmWNWKwYMECm8vrc2zbq1W1R/8MqNGyLjOWYzvq9PJbQs0hrth13333Xa01w2C9TfE94fuy9br65H/Zy62wdY5rSA4Nnv/KK69UqQHQ9z89pw21J/gs+J4tvffee+r1K1asaHAOTUpKijpfoPbTEmotfX19qyzXy4v3dUQOTV3Op7b2EdR84HlLly6ttt+hRs76Sh/v46jvAjmE1jV8qClHzZHle+McjxocPN8e6xoaPBc1Udi3sM6G7rNLKteLGk/Lsr7xxhtquV4jjc+Ez4aaI8tz988//6yehxrB+rD1+27vd6SmHBr44osv1Otwfq6rBnX/eOedd1QUafn366+/SkOhRmTkyJGqzQw1P/of2tHKy8tl6dKl6nlz585V0bitXkEN7eqKLHxEwIgedbgywBUursr++uuvalfMKKcO5QZcQdQEV6HIn8EVDa4uEQGj+yeubt5///1qz2/Tpo26gtchQkXUi4gZV+tGbzdL+P5xlYKrZ+v24IauH58VV3dLliypUjsTFBQkl19+eY2vxRUkrvbt1Shhu1vW8OBqdvDgwWpfANRu4btCLV1UVJT5ebjawNWH/jzA1QeudlDWptClSxd1xff111+bl+G7RS3D5MmT1fbQ9wV0iUT5LPcFXOVgO1huR1twdWl9TNv7qyu9hkV35513qlvL7aeX37LmF7WOOJ7wf9CvEHGFWlpaWqf3Ro0Gai2wjWzVQDpT13jsq6ih0KE2AP9HzdSGDRvM3y9qZVDLZPn96rW9tX2/NUFNF7YXaiUs143zYufOnautG7l2uKJ3hLqcTy33EdSoomxDhgxR/0cNjDXUmlrCOnFOwJW+I74L5A/qtXvYbqdOnVI13tjPLMuD/TY/P79exwzO76hZxRAA1jUjDdlnb7jhhio1kdbbd/369eqzoVbU8tyNXkbY13755Rcxir5f1Gf4hwY1OeEHwJHdd9G8snXrVrtVs3rS7IEDB9QPveWPTGMdOnRIHbTWXXtx8tAft4RmHVsbHVWVdflxQrUcfpDQ1x4naFQBomoYVaQIRHSdOnWqtgPj9YCACCcbI7ebJawfHDmGA36YUeWKIAZVtjhxoLr84osvlrCwsEatG9+3NWxbVC9bfuddu3at9jzsF+iWiBMVqk3x/aFJEGOUIHhAlSqCMTRVOQpO+v/3f/+nmksQiKG7N75bLNdhX0AAgKp1W2pLPMe2xp8jWW9nNH3gOLMcvwjV2gi0V61apZqBLOHzIEhDgIMg9sknn1TV02hyRFCKpkzrRHYdmkzwA+YKY7zg2MS+ZO9Yx483vl80O9R2rDcE1o2La1vHBVg3VWIfrK25rq7qcj5FwIDv/quvvqr2OfWgt67rrK3poi7fBSBZ9ZVXXpHdu3dXCbJxHtchSMA5Bc1z2GZoNkXQWNNQJY4+l8bXsn1rOtchoEETp1FMFT71C+ScojsFfqzwA4bBq2zRdyhnYK93j77x67oOtLfjb+jQoXLOOeeoH27LgMbdtlt9YRvhBwu1V//+97/VDx9qQVDDVRvkcuGKCbVGjQ1+aoMTFK56vvvuO1m0aJEarwTt5bjqxYnMERC4IC8JV+m4csNJEj/0lidG7AsIZrAf2VJbHgdyhGz9ONhimctTH9YnJpy8EazixPnqq6+qoBA/lKjBQeCiD4iI16FGCvkSyLFBQIl8BvygYFlNuV3uAtsC5wtsJ1uw7Rqzbmxj1LLbOr9Zb1/LGpPmOJ/iGENe3D/+8Q9VW4nyoMzY/20NmumIc3RNPvvsM1V7i6AaZcJxh/dE/okekACWo6YX+yu2Lf6QC4kLHgREzcGnibdFU9KDrrrm7DlNQIMrNzTv1PaDjudh50DEXlNtQ30iOiTropYDB4ZlLQ0ib/3xpqTXdKGZwxKSlLHTWX4WJPWBnpDq6O3WUHrSNUZ2RM1SXdX2PeHAx48WfsRwMsCPMppGaoMfSEDVra1kW1yRWsO21ber/p3v2bOn2vOwX+AAs7yKQ80GrsbwhytIJAMjydJeQFPfqmNc9aFWFM1OSD5GsISTqWXtBL4DJDIjabkhPzhYd12bEep6MsR2trxixT6N40zfzvhekYz8448/VrmStNd8gqtj/GHbIpkUyY24akcSqzXsK7gab8hoo43RkGYBBOp6jV9Nx/qWLVtUAFjbe9h73N5yrBvfKb4rR18ENbZpDz9qixcvVjU0lgOt2TqGHaEu3wWCa9TA4ji0/Hy2mvQRoKPZE3/Y93GO+M9//qMSsG2dKy3PpfW9wG0Iy3OddWcVLGvK37/a9g2cv/GbXJ990imGUEUEjipn/OhaQw8FXG0Dqp1x4GHnrukki52xrqNhookAOSmWOQp4P2T440oA1d2OgB5Fttr/9XwC6yo/HFi46teh+hy9Q3CFol8hO3q7NRSqUlETgisUtHHXdf362BX2visEI/hDjwLkAaFXRV3GaEGtl94+bAtyK9B8o0PPBOTB6AEIAhRsZ1xFWZYNJxnUwmCfATQdWtdq4KoM1db4oW7o57ZXS4PaCIxaizZly+YmfV9AeZ5++ulqr8V+UNt7NUUODXLtLOGYAn0761ePlvsItieuYq1/1Kz3I70nib3trE9RgKDJ1n7QVFeo+g9hfb5bfD/4kdOhpxT+j6AMzZj694t91la+HWrX8CNsWQZb72+vbOhVhO8C5wfr7YL/I/+kObeHJVv7CDTVlDl1+S5slQnnD5yLLVlvN+yT+gWWvf0WF0MILPH5rLdZU+yzAwcOVOes9957r0qZcAGJJk7k0jQVvaeTPchZQu4paqObtIYGH1avwbA0bNiwBuUOoNoOV2noqoaqPOw4OEDRLRHRMNoucVWMphl0C0P3T0ToepUjggU8hqtXwOtxtYrqWfy4YAdB0qctyF/BDov3xQZEFI73RBMHdipHNVmgGQLrx8lD36mRQIYgBbUm1vPAICpFV+t169apxGH8kKE7neXJ3tHbraFwJYwmAlwpo+sfmorQVosrSuRF2KteRU0CujojmMTnxXZA27Fl+zFqadDlEOrS3ATYB7EO7AN6V0tLuDIaMWKE6o6JgxjfM5qpLJvu0HSEH14ER/ge9G7bOLj0odzRpIVuuxgSv2/fvioAxnviO0PNkj11+dzW8IOG7YA/PN/66g2BN5IXEVSimhtBJnIf8H2jqQpJ6Chnc+bQ4AoLXY2xv+Fkj6p67BvYVoAy6lewKDtqG/GDjROsZY0l9h80OyJJHlew2O54HvY7Pbi0NwUAAlBsG72rM9aL7YHcgKYY4FP/0XvkkUdUAI7vAJ/POi/DEs5ROD/geMX+gP0C3yG6Juv5Kzh+0dSIhFfUYKEmDgEszsNYjosavbbX3vnPXtmwTZ955hnVrIkyIBDEeQ/fHy6qsO30Y7Ch2wOdLBA0IxiwHGqgNviO0SUYuWq4IEQuCr5TlK0p1OW7wPkWtTPYH/GDj7IgIMAxjX1Yh/MhasVR84HzBPJVcA5BMK7naFpD0IOhDPC94HmoNcVxie8Z3dRtXbw2hp+fn/q8eB8cJ+gco3fbxm+hPs1HU8C+ge173333qd8NnD/xuQHfNTrkoEarXurTJau2rp2W3Wnr020b0MX24Ycf1jp16qQGWMLAb8OGDdNefvnlKsPxo1sbun9169bNPJAZBiCyHMhs9+7darA0vZt0XQbWw0BNeE+sE4MYWXcNrmmobntdcC2hWyU+N7rHoauvn5+fGtgNgz4dOHDA7sB6GFgIXeLxeefMmdOk262h3bZ1P/74o3pvbHd0Vxw0aJDqXmmv+zKsXLlSO+uss1SZbG3HY8eOqS6lXbp00eoDA+dZdmu2/jzomtm+fXu1bTFAmGVXYt3vv/+uhrbXP8/kyZOrDKyH7pAY2AyDH6LbOrrL4r71wHD1+dzW3bYtoSx47G9/+5vdz/3f//5XrRdlRpmwLz/wwAPa0aNHteaifwZsqylTpqhyYKC2O+64o1q3fuwz2Mf1QcxeeOEF88CK+v6F7vroxo/jRR8w8MILL1QDoFmytf8cOnRIdd/G/o7Xojs4jsO6DKzXkG7b+nAKGMgM3axr68JtazA3rPftt9+u9lwcz9g+eD4+C7Ypvusnn3xSDetQl/NfTWXDMBEYIA77Mf5wrsC2wtABtW0Le3DeufPOO9X29/LysjmwnjXr7xFDLKDrNQaZw7kTgzVif7Z+nr0pN+ydrxr6XaD79HPPPacew/eAoTTQzdn6OP/222/VYJDYX3GcY/+9+eab1TmttoH1li9fro0fP958XsExYj10Q02WVK7X+jdD3+7Wv28YgA+fA58Hg2DWd2C9hnTbzsvL066++mr1vVoPrPfrr7+aB0msD6/KQpATQWSMq3X0gvJ0aF7BFQraz9HuXFeoykRNDa7sUMMCuOrC1SpqXxp6xUlERE0LtYTIsbFMu3CZHBoiezCKJarW9dna6wpNQ2hCQvBiqycEERE5H31UZFv5gC7Ry4nIGqaIwFg96NGCaF3vYVAfDz74oPojskUfoLKmXKf6JCQSNTUkKSMvpyYREREO7Vrf3O+J/CK9Q0t9MaAhp4QZVjH2BJIf9d4xRI5UWxI0BkxEDSGRs8A5ER05ajJ79mzVScSV37OhmENDRB4JPYFq6/GCnitEzgJDGOhTMNjTs2dPh/ZYNOI9G4oBDREREbk8JgUTERGRy2NAQ0RERC6PAQ0RERG5PAY0RERE5PIY0BAREZHLY0BDRERELo8BDREREbk8BjRERETk8hjQEBERkctjQENEREQujwENERERuTwGNEREROTyGNAQERGRy2NAQ0RERC6PAQ0RERG5PAY0RERE5PIY0BAREZHL8zW6AETkOUrKKuTHLUdlwfZ0ST6RJwG+PjIsKVpuGpkgrSOCjC4eEbkwL03TNKMLQUTub13KKfnHnC2ScrKg2mNhgb7y5rT+ck63OEPKRkSujwENETW5z1Yfksd/3CHlFZrEhgXIdUM6yFkdWkhWYan8d+lB2XwkS/x9vOXjGwfJ0KRoo4tLRC6IAQ0RNRmcXp5fsFv+89dB9f+L+raRZy/tJWGBfubnlJZXyJ1fbJIFO9JVsLPg7pESHRpgYKmJyBUxKZiImiyY+devZ4KZ+8d3kTem9asSzICfj7e8Pq2fdI4LlczcYnn6550GlZiIXBkDGiJqEm8s3qeak+C5S3vLneM6i5eXl83nBvr5yMtT+woe/n7zUdlw6HQzl5aIXB0DGiJyuP8uPSCv/75P3X/swh5y9eD4Wl/Tt32kTBnQTt1/6w/Ta4mI6ooBDRE51KerD8lz83er+/+Y2FVuHJFQ59feMbaTeHuJ/LknU3YezWnCUhKRu2FAQ0QOM3dDqjz6/XZ1/7YxSXL7OZ3q9foO0SFyfu/W6v5naw41SRmJyD0xoCEih/h12zH5x7db1P3rh3VUtTMNcc3gDur2h01pkldc5tAyEpH7YkBDRI22bF+m3P3VZqnQRK4Y2E7lzdhLAK7NkMQoSYwJkfySclm0I93hZSUi98SAhogaBT2SZn6yQUrKK+SC3q3lX5f1EW8kwjQQAqGL+7VV93/YfNSBJSUid8aAhogabHd6jtz40TopLC2XkZ1j5NUr+4pPI4IZ3UX92qjb5ftPyOn8EgeUlIjcHQMaImqQ4zlFMuPDtZJdWCoD4iPlP9eepSabdISEmBDp1ipMTZWwdF+mQ9ZJRO6NAQ0R1VtRabnM/GS9HM8pViP8zr5+kAT7+zr0PcZ0NU1UiS7cRES1YUBDRPWe0uDBuVtlS2q2RAb7yQczzpaI4KrTGTjCmK6x6vavvZmqpoaIqCYMaIioXj5YnqySdX29veTf0wdIfHRwk7wPZuMOC/CVU/klsjU1q0neg4jcBwMaIqqz/Rl58uLCPer+Y5N7yLCkmCZ7L0xaObKLaf1L2OxERLVgQENEdYJmHwycV1JWIaO7xMq1Q0wD4DUlPY/mrz0ZTf5eROTaGNAQUZ18ufawbDqcpZqBnr+8d4MHzqsPdAWHbWnZkltU2uTvR0SuiwENEdUKUxC8/vtedX/WxK7SOiKoWd4X79M+KkiNQLzxMPNoiMg+BjREVKv3lx6UE3klanyYqwfHN+t7n90xSt2uSz7VrO9LRK6FAQ0R1ehkXrG8v+yguv/AxK4qWbc5DU4wBTRrGdAQUQ0Y0BBRjT5edUgKSsqld9sIOa9Xq2Z/f72GZnNqlhrQj4jIFgY0RGRXQUmZfLIqRd2/dUxSsyQCW0MzV0xogOpdtTU1u9nfn4hcAwMaIrLrq7VHJKugVDpGB8vEns1fOwMIos40O500pAxE5PwY0BCRTRUVmny4Ilnd//uoRIfMot1QAzu2ULfs6URE9jCgISKblu0/IamnCyU80FcuH9DO0LL0ax+pbrccyVJzSRERWWNAQ0Q2fbnmsLq9bEA7CfTzMbQs3VuHi5+Pl5zML1FBFhGRNQY0RFRNRk6R/L7ruLp/1aDmHXfGFgRUPVqHq/ubj7DZiYiqY0BDRNXM2ZAqZRWaDIiPlK6twsQZ9LVodiIissaAhoiqQI7KvI2p6v60s42vndH1bVcZ0KQyoCGi6hjQEFEVO4/lyIHMfPH39ZbzehvTVbumGhpMVFlaXmF0cYjIyTCgIaIqftxyVN2O6xYn4YF+4iwSY0LUTN9FpRWy93iu0cUhIifDgIaIqow989NmU0Bzcb824ky8vb2kT/sIdX/LEY4YTERVMaAhIrP1h07L0ewiCQv0lTFd48TZ9GprCmh2HGVAQ0RVMaAhIrP5246pW0xzYPTYM7b0amMKaLYfzTG6KETkZBjQEJG5uWnB9nR1/3wDZtWuTw3N7mM5UsbEYCKywICGiJStadmSnlMkIf4+MrxTjDijDlHBEhrgK8VlFaonFhGRjgENESkLd5hqZ87pFueUzU16YrA+YjDzaIjIEgMaIlKD6enNTcifMcKOHTtk6tSpEhsbK0FBQdK7d295/fXXpaKiatNSjzamgGZ7Ws15NCtWrJBJkyZJVFSUhIaGyqBBg+STTz6x+dzTp0/Lww8/LOeee6506NBBgoOD1V/Pnj3lgQcekBMnTth83auvviqXXXaZdO7cWSIiIiQgIEC9/rrrrpNt27Y1eFsQUf15aZy6lsjjYVyXCa8tFX8fb9n42HjVrNOcVq1aJePGjZPCwkIVeHTs2FGWLl0q6enpKsj5+uuvxcvLSz332w2pMmvOFhmUECXf3DzU5vrmzp0rV155pQqGRo0aJTExMbJ48WLJysqS+++/X15++eUqz9++fbsKoBD8IIhp06aN5Obmyvr16yUjI0P9f/ny5ZKQkFDldVhvfn6+9OnTR9q2bWsOzPbu3St+fn4yb948ufDCC5tsuxGRBQQ0ROTZ3vh9r9bhwZ+1G2evbfb3Likp0RISEnBhpb366qvm5bm5udrQoUPV8tmzZ5uX7zqWrcra67EFWnl5RbX1nTx5UgsPD1evmzt3rnl5enq61qlTJ7V8yZIlVV6TlZWlrV+/XisvL6+yvLCwULv22mvVay6//PJq77V8+XL1HGvvvPOOek3Lli210tLSBmwVIqovNjkR0ZnmJgN6N3333XeSnJwsffv2lXvvvde8HM1Eb7/9trr/yiuvmJcnxYaqaRlyi8vk8KmCauv73//+Jzk5OXLxxRer5iBdy5Yt5cUXX6y2PkBz0VlnnSXe3lVPiYGBgfLcc8+p+3/88Ue19xo+fLh6jrXbbrtNkpKS5Pjx47Jz5856bQ8iahgGNEQOlpKSoppHxowZo5oj7rvvPmnfvr3KCxkwYID89NNP5ufOmTNHBg8eLCEhIeoH96677lLNLtYKCgrkX//6l/Tv31/90ONvyJAh8vHHH9ssw7Jly+SOO+5QTSEtWrRQ792tWzd56KGHVLOLpdTTBbJxzXI59MKF8t3rj8ipU6fk1ltvldatW6uckF69esmHH34oTeWXX35Rt1OmTKn2GLZXYmKiahLCdgU/H2/pXjkD+A4b49HUtL4LLrhABSC///67FBUV1al8aDoCf3//en2uhr6OiBqGAQ1REykpKVF5IZ9//rkKPvC3ZcsWufTSS9UP6muvvSZXX321hIWFycSJE6W8vFzeeust+dvf/lZlPcjhGDp0qPzf//2fyikZPXq0ygvZvXu3XH/99XLnnXdWe+9//OMf8sEHH6hABmXAH2otXnjhBRkxYoTk5eWZn7tkd4b5fmFejnqvH3/8UUaOHKlqIPA+N910k6r5aArYJnrwYou+fOvWreZlPSvHo9luo6dTTetDcIEADcEM8lxqU1paKk888YQ5GKqrTz/9VPbs2aOShfFHRM2g3o1URFSj5ORklT+Bv7Fjx2p5eXnmx5ALguXI5WjRooW2bt0682NpaWlaXFycevzAgQPm5ZMmTVLL7r77bq2oqKhKTsjAgQPVY7/++muVMsyfP1/lhVjCa2fOnKme/+STT5qXX//hGq3lVc+Zyzxt2rQq7/Pdd9+p5fHx8dU+6+jRo82vq+ufZT4MYDtg+ZYtW2xuz3vuuUc9/uabb5qXfbY6ReXRXPO/1VWem52dbX4f3LflkksuUY//+OOPNh+/8cYbtRkzZmgXXXSR1rZtW/Xc4cOHaydOnNDsefHFF9VrpkyZovXs2VO9pk2bNiovh4iaR/N2ZSDyIMjHePfdd1Vzkg7deVF7sn//fvnnP/8pAwcOND+GnjTTp09XNTfo4YOmls2bN8v8+fPl7LPPVl2ELXM80ET13//+V9VE4H3OO+8882Pnn39+tfKg+QjdoNF89MMPP8hjjz0mhSXlsvLASfNzwsPDVd4Knqu75JJLVK2G3uyDHkg6vKfl/+uiU6dOVf6v1xahm7Qt+vZDryPrKRB2Hs1RXc71HlCWNU/1WZ8lNOOhtkyHpsPZs2dLdHS03c+0cOFC1YtKh67b6CKOvBwiah4MaIiaCH7ou3TpUmUZAhL82GFckwkTJlR7DYIYOHbMNKfSokWLzEGFdcIq6Dk1a9eurfZYWlqaytdBkxGam/TxXNDssm/fPnV/5YETatTdmNAAOS6ifoBt/XDjcyCgQbksAxjk5Biha6sw8fH2kpP5JWp049YRQQ5bd1lZmbrFZ8VYNhifBl26v/32W9U0aAuaEAH5SRh/5qmnnlJNg88884w88sgjDisbEdnHgIaoiejjklhDAGLvcf2x4uJidasnwuJHsaYfRusEV9TmINhADkhNFlfmz/SPbyE7RKRdu3Y2n4c8H8tyORI+Mwa2Q+KzLUistiwDYCTjznGhsjs9V3ak5ZgDGn37AdaHGqe6rM8WJEUjsRi1YwhokK+EmjXLGjdrkZGRKvcItWrIRXr00UdV4Ip1EFHTYkBD1ERs1ajU53HQa1WQyItuwHWxevVqNXgcuiK/8cYbqsmkVatW5mYkNG2h9gFNNXpC8ID4SPmsjmWy9Pzzz6saoPpA0jM+jy4+Pl4FNKmpqapXljUsB9RsWY8YrAKaozlybo+WahkCGHzu7Oxs9boePXrUeX324Hl6kLJmzRoZO3ZsnXo4YWC/DRs2qFoyBjRETY8BDZET02tM0OSEIKWu47rAs88+KzNmzKjyGLqEo6cUIBg4ll0kgX7e0rNyOoH6WrBggfz111/1eg0CLMuABuPPoGfSxo0b1VQF1rAcrIMd5NHM25hWracT1occJLzOOqBBjRWaztB127o5sCYYERgyMzOb9DVE1HDstk3kxMaPH18lSKkL1HaIneYjjHujz3byR2XtzPCkGPH3bdhklH/++adaX33+0HRjSe8OjRwVa5s2bZKDBw+qpGTr5GM9CENicF3X9/PPP6vmOczZZGtAPFuQIIxpD6CutWSgB3r1eQ0RNRwDGiInhkH3ENQgOfX2229Xyb3WULuBmhKdXvOAcWgsc2gwYu2DDz5o/r8e0IztHidGwrg8mCMJnwM9vCxzXfCZwVbt1IM3TZW092+Rgzs3y+n8kipNWmh6Qk8uzKVkOZ4PJpq0tb6vvvrK5mSSGGRw5syZKqhCHo1lryV8J9ju1pNnYptjPCGMRYNxgND0RERNj01ORE7us88+U92j//3vf8sXX3wh/fr1U3kwyBPBYHNHjhyRu+++29xt+4YbblBD+yN3o2vXrip/Az/MqDFA0xV6RB06dEg2HjbV5JzTNU72bk427PMh3wSfEbUmGFUZE1EibwWjHSPXB4m51k1nkJJ8UMpOpYpWVqzyaEZ0NjXxYIJJdE2/4oor1GvRxIWeW+iJhF5IeA8ss4TA5KqrrlK9zBC4oMs3eomh2QpdwZHAbTlBJqCnGLY1mpb03mHovYbACOVGDdBHH32kRokmoqbHgIbIycXFxcnKlSvl/fffVzUJaIbB/zEODX6AMV3CtGnTzM/HD+u6detUbQyCGIz6ixqQp59+WmbNmmVuAkHLU7dWYdImMkhqHzO3aQ0bNkyV+fHHH1fNWKitQTkxZg+CNctAwpYdR7PNAQ1cfvnlKo8G3aaRJI1Rm5FPg+kgbAVHqNVB7yXUuuAPgQ96TKGpa/LkyaqmCMnGltAtG6M3YxsjsEQwgy7xaBpDIIXvxXrMHSJqOl4YXa8J109ETuierzbJ95uPym1jkuSB87qJq3pnyX55aeEeuahvG3nzqv5GF4eIDMQcGiIPg2uYFZWjA1vWargiPTHY1pxORORZGNAQeZh9GXmSmVusumsPiG8hrqxn5RQIySfyJb/YNMIvEXkmBjREHmbF/hPq9uyOUWrEXVcWGxYgLcMDVD7Q7vTqPcCIyHMwoCHyMCv2m5qbhiW5dnOTdS0NejoRkediQEPkQcrKK2TNQVNAM7yT/dmjXUkvPY8mjXk0RJ6MAQ2RB9mali25xWUSHuhrrtlwdT1YQ0NEDGiIPMvKyvyZoUnR4uNd89gurtbTae/xXCkpqzpqLxF5DgY0RB6YPzOik3vkz0C7FkESEeQnpeWaCmqIyDMxoCHyEEWl5bKhcrqDYW4U0GAUYXsTVRKR52BAQ+Qh1qecVk0yrcIDJTEmRNyJHtBgCgQi8kwMaIg8xIoDpvyZYZ2ia50bydX0amtKDN7OGhoij8WAhsjDBtQb7ibjz9iqodl1LEfKKzg9HZEnYkBD5AGyC0plW+U4LcPdKH9GlxATKkF+PlJQUi4pJ/ONLg4RGYABDZEHWHXwpJoeICk2RFpFBIq7QRf07q3D1H0OsEfkmRjQEHmAlZX5M+5YO2OdR7MtlQENkSdiQEPkQfkz7jJ/ky392keq281HsowuChEZgAENkZtLzy6SA5n5goGBhya6x/xNtvStDGiQK1RazhGDiTwNAxoiD2luQpNMRLCfuKuE6BA1R1VxWYXsSeeIwUSehgENkZtbvt/982fA29vLXEuzic1ORB6HAQ2RG9M0TVZWzt/kjuPPWOuv59EcZkBD5GkY0BC5sYMn8iU9p0j8fb1lYMcW4u76xeuJwaY5q4jIczCgIXJjKyubm86KbyGBfj7i7vq2MwU0SILOKSo1ujhE1IwY0BC5sRV6c1Mn9+3dZCk6NEDio4LV/U1sdiLyKAxoiNwU5jTCCMEwzM0Tgi0N7GBqWtuQcsroohBRM2JAQ+SmdhzNluzCUgkL8JU+laPoeoKzE6LU7boU5tEQeRIGNERu3tw0ODFafH0851A/uzL5edOR0xxgj8iDeM5Zjshj52/yjPwZXVJsqLQI9pOi0gpOVEnkQRjQELmh4rJyWVeZQ+LuA+pZ8/LykrM6mJqd1rPZichjMKAhckMbD2WpGorYsADpHBcqnkZvdtKDOiJyfwxoiNy5uSkpWtVYeBo9MXj9odNSUaEZXRwiagYMaIjceP4mT+qubalXmwgJ9veRU/klsuc4J6ok8gQMaIjcTG5RqWxNzfbI/BkdpnoYVFlLs6IyuCMi98aAhsjNrDl4Sg2q1zE6WNpGBomn0ifjZEBD5BkY0BC5mRUHPLu5STessrv62uRTHI+GyAMwoCFyMyv1+Zsqayg8VfdW4RIV4i/5JeWy5QjndSJydwxoiNxIRm6RSoJFx6ahSZ41oJ41b28vGZoYXWXUZCJyXwxoiNzIqgOmH+4erU21E55OT4peti/T6KIQURNjQEPkRpbtM+XPjOjs2c1NulFdTNth4+HTklVQYnRxiKgJMaAhchOapplrIkZ2ijW6OE6hXYtg6dIyVDC23tLKYI+I3BMDGiI3sT8jT47nFEuAr7cMrBz6n0TO6RanbpfszjC6KETUhBjQELlZcxMGlAv08zG6OE5jbFdTQPPnngw1Pg8RuScGNERuNt3BCA8ff8bagA4tJCzQV04XlMpmdt8mclsMaIjcQElZhaw+aOrhxITgqvx8vGVUF1NO0R+7jxtdHCJqIgxoiNwAevEUlJRLTKi/GlCOqprQo6W6XbiDAQ2Ru2JAQ+QGllfmz2DcFQwoR9UTg/18vFTi9P4Mzr5N5I4Y0BC5gWXMn6lReKCfeZA91tIQuScGNEQu7mResWxNNSW7juzM8WfsOa9nK3W7YHu60UUhoibAgIbIxf25J1M0TaRnm3BpFRFodHGc1rk9Wgpa47alZUvq6QKji0NEDsaAhsjF/VE5YNy4ygHkyLaY0AA5u2OUur+IzU5EbocBDZGLd9deutc03cHY7qaePGTfRL3ZaQebnYjcDQMaIhe2PuWU5BaXqe7afdpGGF0cpzexlymgWZdySjJzi40uDhE5EAMaIhe2uLK56ZyuceyuXQdtI4OkT7sIlXP0+y42OxG5EwY0RO6QP9Od+TP1bXZayGYnIrfCgIbIRR3IzJPkE/lqwLgR7K5dZxN7mnKNVuw/ITlFpUYXh4gchAENkYv6Y5epdmZIYrSEBvgaXRyX0SkuTJJiQ6S0XJMllTVcROT6GNAQuSi9yWQsu2vXG5udiNwPAxoiF5SRUyQbDp9W98+r7LlDdadvMwxKWFRabnRxiMgBGNAQuSDULKCnTr/2kdI6Isjo4ric3m0jpE1EoJqhfFnlxJ5E5NoY0BC5oPnbTE0lk3qzdqYhvLy8ZAKbnYjcCgMaIhecjHJN8kl1//xerY0ujsvn0WA8mrLyCqOLQ0SNxICGyMX8tvO4VFRORtk+Ktjo4risszu2kBbBfpJVUCprk08ZXRwiaiQGNEQu5tftpiaS85kM3Ci+Pt5ybuX8V2x2InJ9DGiIXEh2YamsPGBKYj2PzU2NpufRYAoJDVnWROSyGNAQuZDFu46rAeG6tAyVTnGhRhfH5Q3vFC3+Pt6SerpQDmTmG10cImoEBjRELtjcxNoZxwj295XBiVHq/p97OGowkStjQEPkIvKLy+SvvZnqPvNnHGdM1zjzIHtE5LoY0BC50MzaJWUV0jE6WLq1CjO6OG7jnK6miT3RFR5BIxG5JgY0RC7i1+3H1O2k3q3VwHDkGAkxIdIhOljlJmEGbiJyTQxoiFxAQUmZLNmdaQ5oyHEQHI7uYqql4TQIRK6LAQ2RC/hrT6YUlpZL+6ggNaAeOdbwTjHqdkVll3gicj0MaIhcwPzK3k2TerG5qSkMSYwWby+Rg5n5ciy70OjiEFEDMKAhcnJFpeXyx67j6v75bG5qEhFBfmoGblix3zRPFhG5FgY0RE5u6d5MyS8plzYRgdK3nelHl5qu2WklE4OJXBIDGiJXmbuJvZuaLY+G0yAQuR4GNEROrLisXH7faWpumtSbg+k1pbM6tBB/X285nlMsBzLzjC4OEdUTAxoiJ4ZxUXKLy6RleID0b9/C6OK4tUA/HxkQH6nur00+bXRxiKieGNAQObFft1U2N/VqLd7ohkNN6uyOpnmd1qecMrooRFRPDGiInFR5hSaLd5smTJzQs6XRxfGogGbdIQY0RK6GAQ2Rk9p0+LScyi9RXYr1H1pqWv3jI9V4NEdOFUp6dpHRxSGiemBAQ+SkfqscewaTJ/r58FBtDmGBftKjciTmdWx2InIpPEsSOSm9d9O5Pdjc1JwGdqhsdmJAQ+RSGNAQOaGDmXlyIDNf/Hy8ZFTlxInUPAYl6AENezoRuRIGNERO6PfK5ibMMRQe6Gd0cTzKwA6m7vG703Mkp6jU6OIQUR0xoCFyQr/vNPVuOrc7m5uaW1x4oHSIDhYMFrzhEGtpiFwFAxoiJ4OeTesruw2P6x5ndHE8EsejIXI9DGiInMxfezOkQhPp1ipM2rUINro4HunsjqZmp3UcMZjIZTCgIXIyS/eaZnse05W1M0YZWFlDsyU1S0rLK4wuDhHVAQMaIidSUaHJsn2Z6v6oLqbZn6n5JUSHSHigrxSXVcie9Fyji0NEdcCAhsiJ7DyWIyfySiTY38c8Hgo1P8yb1S++hXnEZiJyfgxoiJzI0sramWFJ0eLvy8PTSP3am2be3nQky+iiEFEd8IxJ5ESW7tWbmziYntH6VwY0mw8zoCFyBQxoiJxEfnGZedyTUZ0Z0DhLDc3BE/mSVVBidHGIqBYMaIicxKoDJ6W0XJP4qGDpGBNidHE8XosQf+kYbeo2v5nNTkROjwENkZPlz7B3k/PoX5kYzICGyPkxoCFysvyZ0V04/oyzNTsxoCFyfgxoiJxAWlahpJwsEB9vLxmSyO7azqJ//JmARsPkTkTktBjQEDlJ/gz0bhshYZxd22l0axWuus9nFZSqgJOInBcDGiInCmiGJkUbXRSygGAGQSZwgD0i58aAhshgaMpYfbAyoElkQONsmEdD5BoY0BAZ7MipQpVD4+fjJQMrZ3kmJxwxmAPsETk1BjREBlt10DS7dt92kRLs72t0cchOYvCuYzlSVFpudHGIyA4GNEQGW1mZP4P5m8j5tI0MkpjQACmr0GR7WrbRxSEiOxjQEBmcP6MnBA9hQOOUvLy8qnTfJiLnxICGyECYJygjt1j1phlQOSotOR/OvE3k/BjQEBlIr50ZEB8pgX4+RheH7DDX0DAxmMhpMaAhMtCa5FPqdgi7azu1Pu0ixcvLNKJzRk6R0cUhIhsY0BAZmD+zrjKgGZTA6Q6cWWiAr3RtGabub2QtDZFTYkBDZBBc7afnFImvt5c5R4Ocf+ZtjhhM5JwY0BAZZH2K6YexZ9sIjj/jApDnBBsZ0BA5JQY0RAZZl2Jqbjq7A3s3uYIBld/T1tRsKSmrMLo4RGSFAQ2RwTU0Azsyf8YVJMaESGSwnxSXVcjOYzlGF4eIrDCgITJAdkGp7Dmeq+5z/ibXGWBPHyto4yE2OxE5GwY0RAbYcPiU+aofw+qTa2AeDZHzYkBDZIB15uYm1s64EtbQEDkvBjREBlhfmRDM/BnX0rd9pHh7iRzNLpL0bA6wR+RMGNAQNbOi0nLZcsQ0a/PZDGhcSkiAr3RrFa7us9mJyLkwoCFqZtvTsqWkvEJiQv2lY3Sw0cWhehrQwZRHs4HNTkROhQENkUH5M6idQc8Zci1nVY5HwxoaIufCgIaomTF/xj0Sg3ek5ajmQyJyDgxoiJpRRYUm6yubKs5mDyeXFB8VLNEh/qrZcMdRUy4UERmPAQ1RM9qfmSfZhaUS7O8jPVqbkkvJBQfY05udDnHmbSJnwYCGqBmtTTY1N/WPjxRfHx5+rt7sxMRgIufBMyqREfkzHZg/4w6JwesPnRJN04wuDhExoCEyrocTua6+7SMkwNdbTuSVyP6MPKOLQ0QMaIiaz9GsQknLKhQfby/pVzknELmmAF8f87QVqw6eNLo4RMSAhqj56L2bkAwcGuBrdHGokYYkRKvb1QxoiJwCAxqiZh9/ht213cHQJD2gOaW64xORsRjQEDUT5s+4lz7tIiXIz0dO5ZfIPubREBmOAQ1RM8gpKpXd6Tnq/sDKHjLk2vx9vc/k0Rw4YXRxiDweAxqiZrDx0GlB794O0cESFx5odHHIQYYknml2IiJjMaAhagbrK5ubOP6MmwY0ySeZR0NkMAY0RM1gXWVCMOdvci992kWoaSyyCkplz/Fco4tD5NEY0BA1sZKyCtl8xDTnD2fYdi9+Pt7m73TlAXbfJjISAxqiJrb9aLYUl1VIVIi/JMWGGF0ccrCRnWLU7dK9mUYXhcijMaAharb5m1qomZrJvYzuGmseYK+otNzo4hB5LAY0RE1sbXJlQjDzZ9xS57hQaR0RqGrhOGowkXEY0BA1IfR82XBITwhm/ow7Qq3b6C6mWpq/2OxEZBgGNERN6EBmnpwuKJVAP2/p2SbC6OJQE2FAQ2Q8BjREzTDdQf/2LdTIsuSehnWKUbOoH8zMlyOnCowuDpFH4hmWqDnGn0lgc5M7iwjykwHxkeo+a2mIjMGAhqgJrU3mgHqegs1ORMZiQEPURI5mFUpaVqFqiugfz4DG3Y3uEqduV+4/oQZTJKLmxYCGqImbm3q0DpfQAF+ji0NNrGebcIkJDZD8knJ23yYyAAMaoiaekJLdtT2Dt7eXjO/RUt1fuCPd6OIQeRwGNERNhBNSep4JPU0BzW87j3P2baJmxoCGqAlkW8y+zAkpPcewpGjVvJiRWyybU00TkhJR82BAQ9QENhw+JZomkhgTIrFhAUYXh5pJgK+PjKmc22nRjuNGF4fIozCgIWoCnL/Jc03s2UrdLtqRLhqiWiJqFgxoiJpwhm0mBHse1ND4+3jLwRP5auoLImoeDGiIHKyotFy2pmar+wxoPE9YoJ8M6xSt7i9ksxNRs2FAQ+Rgm49kSUl5hcqd6RAdbHRxyAATepxpdiKi5sGAhsjB9EHVhiRGi5eXl9HFIQOc2yNO8NVvSc2WY9mFRheHyCMwoCFysFUH9ICGzU2eKi4sUAZUTneBMWmIqOkxoCFycP7MpiOm8UeGJpryKMgzTawcZI+jBhM1DwY0RA608fBpNTFhXFiAJMSEGF0ccoI8mtUHT0lWQYnRxSFyewxoiBwIP14wNIn5M56uY0yIdG0ZJuUVmizelWF0cYjcHgMaIgdabc6fYXMTsdmJqDkxoCFykMKSctVlG5g/QzChctTgpfsy1f5BRE2HAQ2RI/NnyiukVXggx58hpWebcGkbGSRFpRUqqCGipsOAhsjB4880Rf7Mrl27ZPr06dK6dWsJCAiQjh07yh133CEnTpyo8XU//fSTjB49WsLDw9XfmDFj5JdffrH53IqKCnnsscekTZs2EhQUpJ67detWm88tKyuT3r17y7Bhwxo0XxG2T23b6KOPPlLPuf76620ut/wLCQlR5UaZH3zwQdmxY0e919sU8D763E5sdiJqWgxoiJx8/Jk//vhDBg4cKF988YVERkbKhRdeqIKad955R/r37y+pqak2X/f666/LRRddJCtXrpThw4fL2LFjZe3ater1b7/9drXnv/DCC/L0009LRESEjB8/XlatWiXnnnuu5ObmVnvuW2+9JTt37lRlMCr5OSkpSWbMmKH+Lr74YunVq5cKZF588UV1/5prrpGcnBxxljwaJAaXllcYXRwi96URUaPlFJZoiQ//onV48Gft8Ml8h603Pz9fa9myJapAtMcee8y8vKKiQps1a5ZaPmHChGqv2717t+bj46MFBARoK1euNC/fs2ePFh0drfn6+mr79u0zLy8pKdEiIyO1vn37akVFRWrZZ599ptb/0ksvVVl3enq6Fh4ert16660N/lxYb22nn9mzZ6vnzJgxo07L9e3y008/aR07dlTPGT16tPpsdX19Uygrr9AGPLVI7RvL92U2y3sSeSLW0BA5qHYG3XMx9kz7KMflz8ybN0+OHz8uXbt2lccff9y8HLUizz33nGp6WrRokWzZsqXK69544w0pLy+XW265RYYOHWpe3qVLF3nkkUdUkxGeo0tJSZGsrCyZNm2aqv2Bq666SgIDA2Xz5s1V1v3AAw+In5+fPPPMM+JssF1QA7VmzRrVBPXXX3/Ju+++a2iZfLy95Nzu7O1E1NQY0BA5wLJ9plyWkZ1jHLreDRs2qNtRo0aJt3fVwxVBBZqS4IcffqjymJ4nM2XKlGrr1Jchv0Z3+vRpdduihWm4fsD7oflJfwzQfPXpp5/Kv/71L4mKct6pHeLi4uSpp55S9998802jiyMTe5kCmkU7jktFRf1zjoiodgxoiBxgWWUPlhGdHBvQ5OfnVws0LEVHm7qHW9bQoKbl8OHD6j5ybKy1b99eYmJi5NChQ+Yck/j4eHW7d+9e8/MQyGRmZpofQ9IwEpHPOussuemmm8TZXXHFFSooO3DggN08o+YyLClGQvx9JD2nSLamZRtaFiJ3xYCGqJGOnCqQlJMFqmkBPZwcKTY2Vt0i+LAlOTm52uN6MIMgCL1/bGnXrl2V17Vq1UoGDBggs2fPluXLl6tg5r777lNBzAUXXKCe895776nmJyQCW9cWOaOwsDBJTExU95HAbKRAPx8Z0y1O3WezE1HTcP6zEpGLNDcNiI+UsEA/h64bTU16E5J1F+20tDT57bff1H3Lnkh5eXnqNjjYfi6PHuhYvu6VV15RNUIjR45UzUno3jxp0iSVk3Ly5El59NFH5cYbb5RBgwaZX1NUVKSCnoay7n5t+XfDDTdIY6EmCiybzYwyoQfzaIiakm+Trp3Ig5qbRnY21aY40oQJE1TNycaNG+X8889XtSM9evSQbdu2yc0336ySe8ERNSYYwwXvgxwZNFsNHjxYrr32WvXYww8/rMabef7559X/Fy9eLHfddZeq+cCYNXgekoyRRFwf6HJtz/79+2XFihWN+kz6GDnOMK/W2G5x4ufjJQcz8yX5RD4nLyVyMAY0RI1QVl4hK/Y3TUKw/kOMnk5o9lm/fr0KMnQtW7aUJ554Qv75z39WybEJDQ1VtwUFBbXm5qBZxlLPnj3NQYsO7/vBBx+o5FrUeKBmaPLkyWqsl7lz56qgBuVArc+rr75ar8+HWqCaHmtsQKPXajlDAjNq7wYlRMmK/Sflj90ZctOIBKOLRORWGNAQNQISPHOKyiQ80Ff6tItskvfo0KGDyl357rvvVC+jwsJCFXhg5GAEO4D/6/QkXjSzIHCxlUejJ8li3bXVcNx+++3Sp08f1QUcUEuEpqZvvvlGdRu/7LLLVG0KlqMrd01NXc0JCc8HDx5U91Gr5QzO6RqnApolDGiIHI45NESN8NeeTHMvFiQFNxVfX1+ZOnWqvPbaayo5984771S1Dghw9OYiHUYT1oOaTZs2VVvXkSNHVM0FghlMh1CTDz/8UNatW6dGFvbx8VHLdu/erWpqEMzokFdTUlKiAhtngYALARnG3sGYNM4AzU6wJvmk5BWbmguJyDEY0BA1wuLdx9Xt2O6mH6rmlJ6eLt9++63quo1aEkt6zyQ8bk1fhmajmiCPBrkzyI/Rx7vRoZbIVhOWs/R+ysjIUPNSwd133y3OIjE2VDpGB0tpuSbLK5PJicgxnOPsQ+SCjmUXyva0HEG+qX7l3RS2b9+umnism4wwfxF6KaF3EhJzLeFHHDUqqM1ZvXq1efm+ffvk2WefVTU+tf3QIzenuLhYzY1kCc1b6EmlD+ZXWloqc+bMUSMMY34lI6FGZv78+SrX6NixY2r+qpkzZ4ozOadyX0GzExE5DnNoiBoIkw1C//aREhNqmi6gKbz88ssqfwa9nTDbNmofMFYMgg10pbbVUwhTJbz00ktqLBl0w8Zkk/7+/mqaBNSuIMG3U6dOdt8TA/UhGMJ7I/nYEnJqMPHllVdeKRMnTlTNTEgMfuihh6oFVk0J20CfMRvNXehajl5aeiIwapaQ14PgzZkg+J29IkWW7MlQAZgz9MAicgfOdaQTuZDFu0zNTeMq5+lpKpdccolqXkKQgV4/6NF03nnnyT333FMld8bavffeq4IWBDbLli1TyzBrN+ZiwtgyNUGOTvfu3dXIwNYwCN/ChQtl1qxZsmDBApWzg/v6VAPNBSMA4w8QSKEcSP4dMmSIXHfddVUSpZ0JejoF+/tIRm6x7DiaI73aRhhdJCK34IUZKo0uBJGrKSgpk35P/SYlZRWy6N5R0qVl1e7PRDWZ+cl6WbTzuNw3vovcNa6z0cUhcgvMoSFq4OjACGbio4Klc5xp3BeiutJzrjAeDRE5BgMaokY1N8UxB4IanBi8JTVLTuYVG10cIrfAgIaonioqNPOV9fgmzp8h99QyPFB6tgkXNPj/WTmWERE1DgMaonracPi0nMgrkbBAXzk7wfgh9cm1m53Q24mIGo8BDVE9/bL1mLod36Ol+PnwEKKGGdM11pyPhTnBiKhxeDYmqmdz06/bTQHNBb1bG10ccmH92reQiCA/yS4sVbk0RNQ4DGiI6tncdDynWMICfGVEE8yuTZ4Dc3+N6mKqpVmym3k0RI3FgIaoIc1NPVtKgK9pskaihhpTGdD8uZd5NESNxYCGqI7Y3ESOptfQYE6wjNyq83URUf0woCGqIzY3kaPFhgVI78qpD5bu5ezbRI3BgIaoAb2b2NxEju7txO7bRI3DgIaoDtCt9ufKgGYSm5vIgcZ0NY1Hs2xvJrtvEzUCAxqiOli2/4ScyCuWqBB/c94DkSP0ax8pkcF+klNUJpuPsPs2UUMxoCGqg+82pqnbi/q2EX9fHjbk2O7bIztX9nbiNAhEDcYzM1EtcotKZeGOdHX/0v5tjS4OuSF23yZqPAY0RLX4dVu6FJdVSFJsiPRpZ+qRQuRIoysTg9l9m6jhGNAQ1WLeplR1e9mAduLl5WV0ccgNxYQGmIPlv9jsRNQgDGiIapB6ukBWHzyl7l/C5iZqlmYnBjREDcGAhqgG8yqTgYckRknbyCCji0NubDS7bxM1CgMaohqmOvh63RF1/8qz2xtdHHJz7L5N1DgMaIhqGHsmLatQwgN95fxeHEyPmmH27cru2xw1mKj+GNAQ2fHV2sPmrtqBfpzqgJpvGgSOR0NUfwxoiGzIzC2W33YeV/enDYo3ujjkIfRRqHccZfdtovpiQENkw7yNqVJWoUnf9pHSvXW40cUhD8Hu20QNx4CGyIqmnUkGvorJwGRU920GNET1woCGyMqa5FNy8ES+hPj7yOS+bYwuDnmYMd0qu2/vY/dtovpgQENkJxkYwUxIgK/RxSEP07ddpLSo7L69id23ieqMAQ2RhayCEpm/3TQRJZOByejZt5fsZvdtorpiQENk4ftNaVJSViHdWoVJX05ESQYZ193U7LSosqcdEdWOAQ2RRTLwV3oy8KB4TkRJhjmnW5z4+XjJ/ow82Z+Ra3RxiFwCAxqiShhufnd6rgT4essl/TgRJRknPNBPhneKUfcXVDaBElHNGNAQVfpqral2ZlLv1hIR7Gd0ccjDnd+rlbr9lQENUZ0woCESkbziMvlp61F1fxrHniEncG73luLtZRo1+MipAqOLQ+T0GNAQichPW45KQUm5JMaGyKCEKKOLQyTRoQEyOCFa3WezE1HtGNCQx0My8GerD5lrZ5gMTM7i/N56s9Mxo4tC5PQY0JDHQzIwqvX9fb1l6llsbiLnMaGHKaDZeDhL0rM5WSVRTRjQkMf7tLJ2ZnKfNtIixN/o4hCZtYoIlLM6tFD3f67M8SIi2xjQkEc7lV8iP281VedfO7SD0cUhquaS/qYhBOZtTDO6KEROjQENebQ564+okYF7t43gyMDklC7s3VoNsrfzWI7sTs8xujhETosBDXmsigpNPltjam66dkgHJgOTU0Iz6NjKGbi/Yy0NkV0MaMhj/bUvU46cKpTwQF81szaRs7q0fzt1+/3mNCmv0IwuDpFTYkBDHuuzVabamakD20uQv4/RxSGy65xusRIZ7CfHc4pl5YETRheHyCkxoCGPhJFX/9iToe5PHxxvdHGIahTg6yMX9mmt7jM5mMg2BjTkkb5Ye1g0TWRk5xhJjA01ujhEtbpsQDvzIHvZhaVGF4fI6TCgIY9TVFouX68zTUR5zRB21SbX0L99pHRrFSZFpRUyd0Oq0cUhcjoMaMjj/Lj5qBp/pk1EoIyr7D1C5OzQC08PwDFVB6bsIKIzGNCQR8GPwIcrktX9GcM6iq8PDwFyrUH2QgN85eCJfFl54KTRxSFyKjybk0dZdeCk7E7PlWB/H5l2NpOBybUgmLlsgGnk4E8re+kRkQkDGvIoHyw31c5MOaudRAT7GV0conrTm51+23VcjmUXGl0cIqfBgIY8RvKJfFm829RV+/phHY0uDlGDdGkZJoMTotQAe5+vPmx0cYicBgMa8hgfVebOIBGYXbXJlekBOWaKzy8uM7o4RE6BAQ15BIzbMaeyq+uNIxKMLg5Ro0zo2UoSYkLUfq0PQUDk6RjQkEf4au1hKSgpV+N4DEuKNro4RI3i4+0lfx+ZaM4LKy2vMLpIRIZjQEMeMZCengyM2hnOqk3uAL2dYkL9JS2rUH7Zeszo4hAZjgENub25G1MlI7dYWkcEyiX9TF1eiVxdoJ+P3DDc1Hz63l8HONAeeTwGNOTWysor1MkeZo5KFH9f7vLkPq4Z3EGNqYSxlRbvMvXgI/JUPLuTW/tl2zE5cqpQokL8OZAeuR2MpYQRr+HV3/ZKRQVrachzMaAht4WT+7+XmGpnbhjWUYL8fYwuEpHD3TwqUcICfGXnsRxZsCPd6OIQGYYBDbmtRTvTZc/xXAnx95HrhnIgPXJPkcH+5qEIUEuDAfeIPBEDGnLb2pnXftun7iNxktMckDu7aWSCRAT5yf6MPPlxS5rRxSEyBAMacks/bzumamfCAn3N43UQuavwQD+5eXSiuZYGQxUQeRoGNOSWPZte/32vuo9ghrUz5CnTIbQKD1RJ8Pq4S0SehAENuZ0fNh+Vg5n5EhnsJzcMZ+4MeYZgf1956Pxu6v47S/bL8Zwio4tE1KwY0JBbQVX7a5W1MzePSpKwQNbOkOe4uF8bGRAfqab5eGHBbqOLQ9SsGNCQTYWFhfLYY49Jly5dJDAwUNq0aSM33nijpKU1LOHw9OnTcvfdd0uHDh0kICBA3d5zzz2SlZVV7bmlpaWyaNEiueOOO6RXr14SHBwsQUFB0r17d5k1a5ZkZmbafI+OHdE121dWPDRODr1wodx2Tic1zYHlX2Ii82nIte3YsUOmTp0qsbGx6rjo3bu3vP7661JRUaH28ccn91TPm7cxTTYcOl3r+lasWCGTJk2SqKgoCQ0NlUGDBsknn3xi9zh++OGH5dxzz1XHMI5N/PXs2VMeeOABOXHihN332bVrl0yfPl1at26tzgE4XnGM1/QaonrRiKwUFhZqQ4YMQd9PrXXr1toVV1yhDRo0SP0/NjZWO3DgQL3Wl5mZqXXq1Em9PjExUa2vZ8+e6v9dunTRTp48WeX5v/32m3oMfx07dtQuu+wybfLkyVpMTIxa1qpVK2337t3V3ue2O+/RIvueq4X0GqeNuXCKNmPGDPMf1oPXXnfddY3ePkRGWblypRYUFKT2ZRyTOJZwPOD/U6dO1SoqKtTzZn2zWevw4M/auFf+1ApLyuyu79tvv9V8fHw0Ly8vbfTo0drll1+uRUZGqvXdf//91Z6/bds29VhUVJQ2cuRI7corr9QmTZqkxcXFqeVt2rTRDh48WO11ixcv1oKDg9VzunXrpo5pHPv4f7t27bQjR444eEuRJ2JAQ9U88sgj6kQzdOhQLTc317z8lVdeUctx4quP6dOnq9fhJFZaWmpefuedd6rlCDisT344Ua9Zs6bK8qysLG3ixInmsln7v3lb1Un8wjeXaeXlphM7lJeXq8AMr0OwROSKSkpKtISEBLUfv/rqq+blOEZxPGD57Nmz1bJTecXawGd+U8fDc7/stLk+XEiEh4er182dO9e8PD093XwBsmTJkmrH4Pr169UxZX0RdO2116rXICiylJ+fr7Vs2VI99thjj5mXI/iaNWuWWj5hwoRGbh0iBjRkpbi4WIuIiFAnmY0bN1Z7vE+fPuoxnNTq4ujRo5q3t7fm7++vTpSWioqKVI0PrhCPHz9ep/WlpaWZa29SUlLMy3cfy9ESHvpZncDXHKxa47No0SL1/LZt21Y7ERO5iq+//lrtx3379q322IYNG9RjvXr1Mi/7bUe6Oh46PvSztj6l6jEBL7zwgnrNxRdfXO2xefPmqccuvPDCOpcPtSx4TYsWLaos//TTT9Xyrl27Vjv+EKTptaebN2+u83sR2cIcGqrWnp6dnS1JSUnSv3//ao9PmTJF3f700091Wt+CBQtU2/7IkSOlZcuWVR5DO/rkyZOlvLxc5s+fX6f1IZcHuQNw9OhR8yB6//x+m2CA1PN6tpJBCVFVXvPZZ5+p26uvvlq8vbnLk2v65ZdfqhyDlgYMGKDyw7Zv3y4pKSlq2bk9WsrlA9oJJuGeNWer5BeX1Xl9F1xwgcqd+/3336WoqG69pfz8TAn4/v7+VZZv2LBB3Y4aNara8YfXDB8+XN3/4Ycf6vQ+RPbw7E5VbNmyxXyCtEVfvnXrVkPWhyRiJCZCq1at1O0364/IupTTatbhRyf3qJbc/N1336n711xzTZ3eg8gZNeRYemxyD2kdESjJJ/LloXnbUCNfp/UhKEFCPoKZvXtNvQZrgkT+J554whwMWcrPz1e3LVq0sPna6OjoKuUhaijfBr+S3NLhw4fVbbt27Ww+ri8/dOiQIet75513pKysTPXsSEhIkMzcYnlu/i712H3ju0jbyKAqz//+++8lNzdX+vTpo/6IXFVDjiVMh/DWVf1l2n9Xy09bjkr/9pFq3qecnBxVE1vb+tavX6/WZ+vYuemmm1TtKi4wUAuDHpCobXnxxRerPE+vUbV3jCcnJ9frHEBkD2toqIq8vDx1i66YtoSEhKhbBAnNvb5NmzbJM888o+6/8MIL6vbZX3ZKTlGZ9GwTrkZKtfbpp5+q22uvvbZO5SVyVg09lgZ2jJKHJ3VX95/5Zacs3nXcvK6GrE/38ccfq78ff/xRBTNjxoxRzbt6jYsOTU16E5d1F2287rfffqvxfYjqigENuYTjx4/LZZddpqrAMX7N+eefLwu2p8v3m4+Kt5fIvy7rLb4+VXfnjIwMdbJEuz3yZ4g81Y3DO8qUs9qpPLPbv9go29Oqj/9UX6gpRRMWctnmzJkjqampquZ04cKFVZ43YcIE1ayFIArH7dq1a9X9VatWqf9jPcD8Nmos7kFUBQbWgoKCApuP6+3hYWFhzbY+XLlh4C8kO2JAsVdeeUWOZRfKQ/NMuQIzRyVJn3aR1V731VdfqZPluHHjVDIxkStrzLGEAfcQ9I/uEitFpRUy67s95scae6xjoDwkFiOBGO9z/fXXm1+rv/e8efPU4Htowho8eLBa57Bhw9RFh557Yy/HhqiuGNBQFfHx8eoWV1u26MsxSmhzrA81MhdddJFs3LhRXemhSlsTL7nv6y2SVVAqvdtGqNwZW/TeTUwGJnfQ2GPJz8db/j19gPRpFyE5Ff7iExji0GMdz0NvxvT0dFmzZk21xzZv3izffPONqmG9+eab5c0335SdO3eaez8i4CFqDCYFUxV9+/ZVtwggbNGX1zXBtjHrQ+3KlVdeKX/++ae6msNVHnpfYOK9VQdPSpCfj7wxrZ/4+1aPy9EzY926dSo/AE1VRK4OxxJ6AuGYQY1lQ47NkABfmX392TL9f2skPTZByo9sl5//WCE9evSo1msJXcDRdRvTn9RVTEyMurU1PYmvr6+qYcWfpZUrV6pb5OAQNYrN0WnIY1kOrLdp0yaHDqxnPXheTQPrYRTRa665Rr1Xv379tNOnT6vli3akq4HCMGDYV2sP2X3fRx99VL0WoxQTufvAehgE03pgvZpgJOGuF92iXhPebai2PuVUowfWKysrU1Ob4HXr1q2r02uOHTumhYWFadHR0VpBQUGd34vIFgY0ZHfqg2HDhml5eXl1mvrgrbfeUiOBPvTQQ3anPsCQ6JZTH9x11102pz6wfAzzvmRkZKhlO9Kyte6P/qqCmYfnbTXPW2OLfmJdsGBBg7YBkatMfYBj1HrqA0tjx45Vx6b1VCLJR9I138AQ05xtlz+izd96VC3HxYW9qQ++/PJLbevWrTanUbjxxhvVa3r37l3t2MQcUJgewXpkYX2OuI8++qiBW4XoDAY0VA1OPIMHD64yOaX+f3uTUz7++ON2gxNMTpmUlKQexy0mtMOVJP7fuXPnapNTfv/99+bpDcaPH6/WecXV12ixAyaoiScThk7Srr3uOm3Xrl02y79ixQrzJJa4aiRyF9i39ckpcUzi2NTnKZsyZYrNIL9Dhw42gxP47MuvNS8vb03ESwuM76MNGH2+eXLK++67r9rzcSzqk8xiyoSrrrpKGzVqlBYaGmqeXmTnzp02X4d5o8aMGaNeM27cOC0gIEC9BrWpRI7AgIZsQvUvTjQIQNBchODg+uuvtzsrbk0BDSBowWSU7du3V+vDLWph9KYkS7jK1AOamv5snaDh1ltvVY/fe++9jdwKRM5n+/btqrYTzTSBgYFq5nrU2Nibp6ymgAb+/GupltR/uOYdEKJ5+QVorZJ6au9/8KHN5y5btky77bbbVLNXTEyM5uvrqwKgIUOGaM8++6yavNKW7777Tk0si/OIn5+fmp0bAZG9MhE1hBf+aVwWDlHTOZpVKFe9v1oOnSxQowB/NXOItI+yPRAYETUMfgY+WJ4sz87fpeZ+6h8fKe9OP0taRQQaXTSiOmNAQ05r59Ecmfnpekk9XSjto4Lky78PkXYtGMwQNZUluzPk7q82qdG3Y0L95eWpfWVM1ziji0VUJwxoyCn9sDlNHpy7VQ0C1jE6WD7/+5Bq8zQRkeMdOpkvN3+6QXanm6YimD44Xv5vUnfV5ZvImTGgIadSUFIm/5q/Wz5dbZqoblSXWHlzWj+JDPY3umhEHqOwpFxeXLhbZq9IUf9v1yJI/nlBd5nYs5Ua+ZfIGTGgIaexZE+GPP7DDjl8yjQU+21jkuT+CV3FB5M1EVGzW7n/hMyas0WOZhep/w+Ij5Rbx3SScd3ixJvHJTkZBjRkuH3Hc+WZX3bJX3tNo4u2iQiUF6b0kZGdY40uGpHHQ63pe38ekP8sPSjFZRVqWee4ULl6cLxc3K+tRIWw9pScAwMaMgR2uzXJp+R/y5Jl8e7jqmeFn4+XXD+so9w5rrOEB/oZXUQispCRUyQfrkiRz1cfktxi0wzZOGbHdWspl/RvK+d0i5UAXx+ji0kejAENNavS8gqZv+2YCmS2pWWbl5/Xs5U8dH436RhjmjCPiJxTTlGpfLcxTeZsOCLb03LMy8MDfeWCPm3kkn5t5OyOUWySombHgIaaRXZBqXy57rB8tCJF0nNM7fGBft4y5ax2csPwBEmKDTW6iETUgKEVvt+cpnolHs8pNi9Hj8RrhnSQqwa1Z0I/NRsGNNSkDp8skA9XJMs3649IQUm5WhYbFiAzhnaQqwd3YPs7kRsor9BkzcGT8t2mNFmwPd3cJIWLlssGtJObRvCihZoeAxpqEqmnC+Stxfvl242p6mQH3VqFyd9GJsrkvq3Z1k7kpopKy+XnrcfUyMO7jpmapND6dGn/dnL3uM4SH83BMalpMKAhhzqWXShv/7Ff1ciUlmvmsWRuHpUow5KiOYYFkYcl/r+/9KAs3p2hlvl6e8nUge3lzrGdpA0HyiQHY0BDDusB8e8/D8gXaw5LSbmpa+eITjFy7/jOclaHKKOLR0QG2nwkS179ba8srRyawd/HW3X7vu2cJIkL43xR5BgMaKhRTuQVqzEqMLKvPkbFoIQouW98FxmSGG108YjIiaxNPiUvL9qjbiHIz0dmDOsot4xOZPIwNRoDGmqQU/kl8t+lB+XjlSlSWFpuHkUUI/u6etNSZqbpKpKoOcXGesZAkvjJWbH/pApsUHMDYQG+Kr/uhhEdOQYVNRgDGqp309L7yw7KZ6sPmwOZvu0i5N7xXWR0l1iXDmR07vAZyPV42qkYn3fxrgx55be95uThEH8fueLs9nLDsAQmD1O9MaChOtl7PFc+WZUi36xPlZLKpqVebcPlnnFdZFz3OLcKAtzps5Dr8NRTcUWFJr9uT5c3F++TPcdNM3zD4IQoNQIxBt1sweEdqA4Y0FCNg+Et2pkuc9anytoUU5u33rSE6QnGuEmNjDV3/Ezk/Dz9VIzPv2zfCfnf8mRz8jDgcOzZJlyGJ8VI3/aR0qN1uMRHBXMkYqqGAQ1VcTSrUJbvPyG/bjumbvWu15jx+tzucSqBb2iia+fI1MadPxs5L56Kz0jLKpQfNx+VH7ccNTdHWULTVJdWYWqSzE6Vf53jwtQIxQx0PBcDGg+Gqt6DJ/Jla2qWrEs5LasOnJCUkwVVnoPB8C7o3VqNHdEqwjO6VzIpmIzgKUnBDcnbW3ngpKw+eFJ2HsuR3em55mZvaxiZGCMSqyAnNlQ6tzTd7xAdIn4+3s1edmpeDGg8aPTO/Rl5KhdmT3qubE3Nlu1p2eYhynW4uOndNkLGdW8pk3q3VicDIiJnUVZeoS7EcC7bdzxP9mfmyf7jeZJ8It88BpY1DOiHiW8R5OCclhgbIomxplv2qnIfDGjcTG5RqRw+VSCHThaoAx5/uKJJOZEvlTMQVLui6dUmQvq1j5ShSdFydkIUD3AicslA58jpQtl3PNcc5KjbjDzzPHK2xIQGqMAmCUFOTKjqXdUqPFBahgdKTKi/+LJmx2W4dUCTX1wm8zalib+Pl6pu9Pf1Nt36mG79fLzUzopaCeSIeHuZ/nC/QtNUt+SiknJ1q/5KytWBkV9SptadX4z/m26xrKC4XPKKy6S0vEKtw9cH6/JWVwcBvt4S6OdjvkUggfmMAvy8JdDi1vo5KDPmQkIuCw5YXIHgfnZhqZzKL1bjwZzMK5GM3GIVyOD/9kQE+UnXVmHStWWY6qHUp12kaoPmAUtE7ty0fiynSAU2+t/BzDxVy5OZe2aGcFuQThcdEiCRwX4SFuirLvZwGxboJ6EBPmpgwCB/Xwny85Ygf/2+jwT7m87l+n3TY6b/O6LpCz/bZRWa+m3Ab0JxaYWqhcfgpsVlplvz/0srl5VWqN+Zi/u1FXfl1gENJkgc8cIS8TSYwRq9AFC1ihyYLi3D1C1muWbCKxHRmRptNFUdzMRfnhw4kS+ppwtV3g4uEvWJdR0JF7gIbnDRi7MxzsmmWzzqpW7P/N80kzkuYk23FeZApiHaRwXJsgfGirvyFTeG2g2MYVBqrtkw1W6o/5eZ/so1086BsA63qJnBH3apIH9vUwSuaktMf6EBviriDgnwlZAARN++5mXqNsBX1fxUVIiUVVSYd0Y9Qi7Sb1H7U1YuRZXRM271iFrdVj4XZcSOb65Z8vUSX29vdZUQHeIvUSEBEhXqL7Gh/tI+KlgFMrh6ICKimuFciZpq/FnDuftkfrFk5BRLTlGp5BaVSU5h5W1RqbnGXq+9N9fiV9bsF5SWSWFJhRSWlKnH9BgEAQnW4WiBeq2/RWtAQOUy/bG4sABxZ25dQ0NERGQ0/MzioloPfBAIoSkMP774BcY9dWtxX2e6oDVdyOoXt0hn8MP/VTqFl7rY9WLtOwMaIiIicn3MBiUiIiKXx4CGiIiIXB4DGiIiInJ5DGiIiIjI5TGgISIiIpfHgIaIiIhcHgMaIiIicnkMaIiIiMjlMaAhIiIil8eAhoiIiFweAxoiIiJyeQxoiIiIyOUxoCEiIiKX5yvOMrV6SYnRxSAiIqIG8Pf3Fy8vLxFPD2gQzDz//PNGF4OIiIga4KGHHpKAgAAxkpeG6hEnrKFJT0+Xjz76SK6//npp1aqVeDJuizO4Lc7gtjiD2+IMboszuC2ab1uwhqYSNoJ1ZIeNo98aHfUZjdviDG6LM7gtzuC2OIPb4gxuC8/aFkwKJiIiIpfntAFNaGiojB49Wt16Om6LM7gtzuC2OIPb4gxuizO4LTxrWzhFDg0RERGRW9bQEBEREdUVAxoiIiJyeQxoiIiIyOUxoCEiIiKXx4CGiIiIPDegWbdunUyaNEkiIyMlJCREhgwZIt9880291lFcXCxPPfWUdO7cWQIDA6VNmzYyc+ZMycjIqPbcgoICeeWVV+Tqq6+Wbt26ibe3txqQLyUlxe76x4wZo55j669jx47iKM29LTZv3iyPPvqoep+4uDg1SFJiYqLcdtttkpaWZvc99u7dK1dccYXExMRIUFCQ9O3bV9599101UrMnbQt33S/27Nkjf//736V///4SGxurtgU+z4UXXiiLFy/2qP2iIdvCXfcLW/D++Fx4rSftFw3ZFu66X6SkpNj9XPh74oknbL7HsWPH5KabbpLWrVur9+jatas8++yzUlpaKkZr0EjBS5YskYkTJ6oPM23aNAkLC5O5c+fKlVdeKUeOHJH777+/1nVUVFTIxRdfLAsXLlRf3OWXXy779u2T//3vf+qEs3r1anUi0uELmTVrlrrfoUMHadGihZw6dapO5X388cerLcNO4whGbItbbrlF1qxZI4MGDVLviZM1/o+TzZw5c2TZsmUq6LO0c+dOGTZsmBQWFqqTFHb0X375Rf3w47G33nrLY7aFu+4X27Ztk3nz5snQoUPVdx0eHq6Cuh9++EF9188884w88sgjHrFfNGRbuOt+Ye39999Xr0UZ7AUn7rpfNGRbuPt+0bdvX7nkkktsBnK2pk8YPHiwpKamyqWXXqoCp7/++kv++c9/ytq1a+X77783dvoDrZ5KS0u1pKQkLSAgQNu0aZN5eVZWltalSxfN399fS0lJqXU9H374IfYe7aqrrtIqKirMy9999121fObMmVWen5ubqy1atEg7efKk+v/EiRPV85KTk+2+x+jRo9VzmopR2+LNN9/U9u3bV209zz//vHr+pEmTqj02atQo9dj8+fPNy4qLi7WRI0eq5StXrtQ8ZVu4635RVFRU5Xm6tLQ0LS4uTvPz89NOnz7tEftFQ7aFu+4XlnC+DAsL02bNmqV16NBBlcUWd90vGrIt3HW/SE5OVstnzJhR57Jed9116jVYpw7vNW3aNLX8iy++0IxU729p4cKFquA33HBDtcc++ugj9diTTz5Z63qGDh2qnmv9RWHjJCYmaiEhIVpBQYHd1ztDQOMs20JXVlamBQUFqedb2rNnj1r/OeecU+01f/75p93P4I7bwhP3C7j00kvVujZv3uzR+4W9beEJ+wUex3eNH0k8bu9H3BP2i7puC3feL5LrGdDk5OSobYR1WV8s4D3t7TPNqd45NH/++ae6nTBhQrXHUGUGqIKqSVFRkWoWQNsbmo8sobpq/Pjxkp+fL+vXrxdH+OKLL+S5556T119/XZUfVXOO4GzbAs/38/MTX1/fOpdzxIgRqr22tnK6y7bwxP3i5MmTal3BwcEqv8iT9wt728IT9gs0E2H9H374ocqJaUg53WW/qOu28IT94ujRo/LOO++oz/bBBx/IgQMHbL7HqlWrVI4O1mXdrIT3xHuvWLFCysvLxWVyaNAmB2g7s4YpyTFPhP4ce7DBsDPYWoflurGekSNHSmNNnz69yv+7dOkin3/+uQwcOLBR63W2bfHtt99KTk6OTJ06tc7l9PHxkYSEBNUuXlZWVmMA4A7bwhP2CyRz4uSLEwtOVj/++KNkZWXJe++9p9rmPWm/qOu2cPf9Av9/+OGH5a677pLhw4c3uJzusF/UZ1u4+34Bv/32m/rTIVjBZ8UxguC1LuXUlyMZ/9ChQ3YvFppavWtosrOz1W1ERITNx5F8pz+nMeuwfF5DIUHq559/VsmA6CWFg/Duu+9WXz6izMOHDzdq/c60LZA4hgMUVxtPP/10vd8DB0Nubq64+7bwhP0CP+JPPvmkSnzFFSiu3mbPnq16JnjaflHXbeHO+wW+wxkzZqheKeiN4ohyuup+Ud9t4c77RXBwsOohumHDBhXko5PN77//rjpYfPbZZ3Ldddc1+j2aW8PCaxdx7733Vvl/9+7dVXUhNjx+6F5++WV58803xdWhGh3d/dAT7JNPPlFVf56qLtvC3fcLdE1GflxJSYnqmomeHDg5oReCK3+upt4W7rpfvPTSS6qHC3rS4EfMkzVkW7jrfhEXF6e6eFsaN26c6hk4YMAA1Utw48aN6r6rqHcNjR6d2YvCUM1vL4Krzzosn+doN998s7pFe19jOMO2wA84dsIdO3aorsrXXHNNg94D1Yz2quDdaVt4yn4B/v7+qlocJ/Fbb71V5Q38+uuvHrdf1GVbuOt+gRoqdDdGd+vRo0c7rJyuuF80ZFu4635REwR61157bbXPZvTvdpMENJZtcbb6qOfl5dltY9OhfQ0D49lrF6ytra6xoqOj1QGJJKnGMHpb6D/gW7Zskbffftt8gNWnnMgtSE5OVu3iDW0Pd6Vt4Qn7hS16wqGegOgp+0Vdt4W77hdoHkEiJ5I+rQdOQ64DHtP/j2YHd94vGrIt3HW/qA0GUwTLz1ZTOfXluHCIj48Xlwlo9Mh20aJF1R7DgD6Wz7EHuQ1op9MTiCyhehgJSkhGamyylT2obsb7NHaURyO3heUPOK42cdXRkHIuX75c7bSNvWJxlW3h7vuFPUiIBfT88pT9oj7bwl33C5QZ+UK2/pBsiiRf/f8YlNKd94uGbAt33S9qgx5TYPnZMFgfAhasy3oQQrwn3htJ1o0JdButIYMAoR96TYMAWY4Nc/ToUW3Xrl3qcUcNiFSXcWgOHjxoHoTPUmpqqtazZ0/12o8//lhrDKO2BT5Xv3791GNvvPFGncpa20BZK1as0DxhW7jzfrF+/Xqbg8lhjIj27dur1yxfvtwj9ov6bgt33i/saczAeq66X9R3W7jzfrFx40abx8jcuXM1b29vrUWLFtXew97AenhPlxxYD/744w810iZGWfz73/+u3XfffWqHwAd6+eWXqzwXg/Zg+ezZs6ssLy8vNwclQ4YM0R588EHt8ssv17y8vLSEhAQtIyOj2vvef//9an34a9OmjXotXqMvW7Zsmfm5eL/AwEBt/Pjx6ovE+q+44go1uBBeN336dJtfpitsC32gp27dummPP/64zT/rUVC3b9+uRUREqIPj2muv1R544AHzAXnHHXc0eju4yrZw9/2ibdu22mWXXabdc8896j0vvvhi9Z1jHRgR1Zo77xf12RbuvF80JKBx1/2ivtvCnfeL0aNHa+3atdOmTp2q3Xvvvdpdd92ljRgxQr0e2+KHH36oVk4EU7ggwDqxbrwH3guvmTx5skO2RWM0ePjDNWvWaOedd54WHh6uRmQdNGiQ9tVXX1V7nr2Nrw9P/sQTT6hhn3HgtGrVSvvb3/6mpaen23xP/Qu292f5Hlu2bFEHYo8ePbTIyEjN19dXi4mJ0SZMmGCznI3R3Nuitu1gr+Zq9+7d2pQpU7SoqCi1w/bu3Vt75513HLoTOvu2cOf9Ys6cOerkpI8KihOk/qO+YMECu+V0x/2ivtvCnfeLhgQ07rpf1HdbuPN+8f7776v3Q4CC99NHAcbzUQNkD4KaG2+8UWvZsqV6j86dO2tPP/20qsEzmhf+Ma7Bi4iIiMiApGAiIiIiZ8OAhoiIiFweAxoiIiJyeQxoiIiIyOUxoCEiIiKXx4CGiIiIXB4DGiIiInJ5DGiIiIjI5TGgISIiIpfHgIaIiIhcHgMaIiIicnkMaIiIiEhc3f8DgFroQt0nV4AAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "az.plot_posterior(trace, var_names=[\"beta1\"], hdi_prob=0.90)\n",
    "plt.title(f\"Effect on {metric} (positive = {model_2} better than {model_1})\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "debcca8b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(model_2 > model_1) = 1.000\n"
     ]
    }
   ],
   "source": [
    "beta1_samples = trace.posterior[\"beta1\"].values.flatten()\n",
    "prob_improvement = (beta1_samples > 0).mean()\n",
    "print(f\"P(model_2 > model_1) = {prob_improvement:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fc5007a",
   "metadata": {},
   "source": [
    "# Sequential Bayesian Hypothesis Testing with Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "d16663f5",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d8c8815d96314e3ca504e83217cdb9cd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 39 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(a)] Accuracy: cached -> bayes_cache/post_a__Accuracy.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7fefb15749f0402e91771d38f28390ad",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 26 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(a)] Precision: cached -> bayes_cache/post_a__Precision.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a61938b6f8c4410db44496b78e179680",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 26 seconds.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(a)] Sensitivity: cached -> bayes_cache/post_a__Sensitivity.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5f418a8df4d642f887608171a37dd1c6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 29 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(a)] Specificity: cached -> bayes_cache/post_a__Specificity.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4dcd760478cb4e29a8f36055b56595a3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 27 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(a)] F1 Score: cached -> bayes_cache/post_a__F1_Score.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "37c7095bcecf404b97f9c842248cb01c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 28 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(a)] MCC: cached -> bayes_cache/post_a__MCC.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bec8e84f876a4e598cb59c8e0d05c4d6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 29 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(b)] Accuracy: cached -> bayes_cache/post_b__Accuracy.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0e8be6c079364285aad6c1fcca78386c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 27 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(b)] Precision: cached -> bayes_cache/post_b__Precision.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "018b8115ad0a4d2e931a066289d1a16c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 26 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(b)] Sensitivity: cached -> bayes_cache/post_b__Sensitivity.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a319c073886f41a59947115ec63abc12",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 25 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(b)] Specificity: cached -> bayes_cache/post_b__Specificity.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bed969a076e34a6da707979d63004882",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 25 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(b)] F1 Score: cached -> bayes_cache/post_b__F1_Score.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2d5ed31fa4e34853a1778539e8fff7fd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 24 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(b)] MCC: cached -> bayes_cache/post_b__MCC.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f420bca730e14b64a628d94cfd53ff59",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 25 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(c1)] Accuracy: cached -> bayes_cache/post_c1__Accuracy.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7930cb5228144d758a4d4eed88b100ac",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 27 seconds.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(c1)] Precision: cached -> bayes_cache/post_c1__Precision.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a2b8c069c3664e4eb25605092a047645",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 27 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(c1)] Sensitivity: cached -> bayes_cache/post_c1__Sensitivity.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2df3661fce7d4ac288fa855c3db4ddd6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 26 seconds.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(c1)] Specificity: cached -> bayes_cache/post_c1__Specificity.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ed3dcc8af5fa405cb6a4c415565abc8a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 25 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(c1)] F1 Score: cached -> bayes_cache/post_c1__F1_Score.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "55af1e8ea8ac41e8a065635a5340a88c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 26 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(c1)] MCC: cached -> bayes_cache/post_c1__MCC.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0a9ba2c85ee5487784a8361cdf01b28f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 27 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(c2)] Accuracy: cached -> bayes_cache/post_c2__Accuracy.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1bffaf97b60246f09982b5cadf45312c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 26 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(c2)] Precision: cached -> bayes_cache/post_c2__Precision.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "65dc9080e0914f8b9222738ee05b9d6f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 25 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(c2)] Sensitivity: cached -> bayes_cache/post_c2__Sensitivity.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fa1bf4ffaede45b78f40d132c24c583c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 25 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(c2)] Specificity: cached -> bayes_cache/post_c2__Specificity.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "00e0ce26547f4ab9aa0a5fbbc5ff9f95",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 25 seconds.\n",
      "Initializing NUTS using jitter+adapt_diag...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(c2)] F1 Score: cached -> bayes_cache/post_c2__F1_Score.npz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [beta0, beta1, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6b8d662dcc224a528662c13cfc900a92",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 26 seconds.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(c2)] MCC: cached -> bayes_cache/post_c2__MCC.npz\n",
      "Done. Cache written to bayes_cache\n"
     ]
    }
   ],
   "source": [
    "# ---------- config ----------\n",
    "CSV_PATH = \"results_metrics.csv\"\n",
    "CACHE_DIR = \"bayes_cache\"\n",
    "os.makedirs(CACHE_DIR, exist_ok=True)\n",
    "\n",
    "metrics = [\"Accuracy\", \"Precision\", \"Sensitivity\", \"Specificity\", \"F1 Score\", \"MCC\"]\n",
    "\n",
    "# (label, MODEL_1, MODEL_2)  -> Model_1 is baseline (group=0), Model_2 is treatment (group=1)\n",
    "comparisons = [\n",
    "    (\"(a)\", \"basic_no_fft\", \"cdd_fft\"),\n",
    "    (\"(b)\", \"basic_no_fft\", \"cdd_no_fft\"),\n",
    "    (\"(c1)\",\"basic_no_fft\", \"basic_fft\"),\n",
    "    (\"(c2)\",\"cdd_no_fft\", \"cdd_fft\"),\n",
    "]\n",
    "\n",
    "# MCMC settings\n",
    "draws = 2000\n",
    "tune = 1000\n",
    "chains = 4\n",
    "target_accept = 0.95\n",
    "random_seed = 2025\n",
    "\n",
    "# ---------- data ----------\n",
    "df = pd.read_csv(CSV_PATH)\n",
    "df = df[df[\"Dataset\"] == \"External\"].copy()\n",
    "df[\"Model\"] = df[\"Model\"].replace({\n",
    "    \"Basic Model No FFT\": \"basic_no_fft\",\n",
    "    \"Basic Model With FFT\": \"basic_fft\",\n",
    "    \"CDD-based Model No FFT\": \"cdd_no_fft\",\n",
    "    \"CDD-based Model With FFT\": \"cdd_fft\",\n",
    "})\n",
    "\n",
    "# ---------- sampling ----------\n",
    "index_rows = []\n",
    "for label, model_1, model_2 in comparisons:\n",
    "    lab = label.strip(\"()\")\n",
    "    for metric in metrics:\n",
    "        y1 = df[df[\"Model\"] == model_1][metric].dropna().values\n",
    "        y2 = df[df[\"Model\"] == model_2][metric].dropna().values\n",
    "        if len(y1) < 2 or len(y2) < 2:\n",
    "            print(f\"[{label}] {metric}: skipped (insufficient data)\")\n",
    "            continue\n",
    "\n",
    "        y = np.concatenate([y1, y2])\n",
    "        group = np.concatenate([np.zeros_like(y1), np.ones_like(y2)])\n",
    "\n",
    "        with pm.Model() as m:\n",
    "            beta0 = pm.Normal(\"beta0\", mu=0.5, sigma=1.0)\n",
    "            beta1 = pm.Normal(\"beta1\", mu=0.0, sigma=1.0)\n",
    "            sigma = pm.Exponential(\"sigma\", lam=10.0)\n",
    "            mu = beta0 + beta1 * group\n",
    "            y_obs = pm.Normal(\"y_obs\", mu=mu, sigma=sigma, observed=y)\n",
    "\n",
    "            idata = pm.sample(\n",
    "                draws=draws, tune=tune, chains=chains, target_accept=target_accept,\n",
    "                random_seed=random_seed, return_inferencedata=True, progressbar=True\n",
    "            )\n",
    "\n",
    "        b0 = idata.posterior[\"beta0\"].values.reshape(-1)\n",
    "        b1 = idata.posterior[\"beta1\"].values.reshape(-1)\n",
    "        p_gt0 = float((b1 > 0).mean())\n",
    "        hdi95 = az.hdi(b1, hdi_prob=0.95)\n",
    "        hdi90 = az.hdi(b1, hdi_prob=0.90)\n",
    "        #hdi90_low, hdi90_high = az.hdi(b1, hdi_prob=0.90)\n",
    "        #hdi95_low, hdi95_high = az.hdi(b1, hdi_prob=0.95)\n",
    "\n",
    "        summ = az.summary(\n",
    "            idata, var_names=[\"beta0\",\"beta1\"], kind=\"diagnostics\", round_to=None\n",
    "        )\n",
    "        # rows indexed by param names\n",
    "        rhat_b0      = float(summ.loc[\"beta0\",\"r_hat\"])\n",
    "        rhat_b1      = float(summ.loc[\"beta1\",\"r_hat\"])\n",
    "        ess_bulk_b0  = float(summ.loc[\"beta0\",\"ess_bulk\"])\n",
    "        ess_bulk_b1  = float(summ.loc[\"beta1\",\"ess_bulk\"])\n",
    "        \n",
    "        \n",
    "        # cache file: one npz per (comparison, metric)\n",
    "        fname = f\"{CACHE_DIR}/post_{lab}__{metric.replace(' ', '_')}.npz\"\n",
    "        np.savez_compressed(\n",
    "            fname,\n",
    "            beta0=b0, beta1=b1, p_gt0=p_gt0,\n",
    "            hdi95_low=float(hdi95[0]), hdi95_high=float(hdi95[1]),\n",
    "            hdi90_low=float(hdi90[0]), hdi90_high=float(hdi90[1]),\n",
    "            rhat_b0=rhat_b0, rhat_b1=rhat_b1,\n",
    "            ess_bulk_b0=ess_bulk_b0, ess_bulk_b1=ess_bulk_b1,\n",
    "            label=label, model_1=model_1, model_2=model_2, metric=metric\n",
    "        )\n",
    "        index_rows.append({\n",
    "            \"file\": os.path.basename(fname),\n",
    "            \"label\": label, \"metric\": metric,\n",
    "            \"model_1\": model_1, \"model_2\": model_2,\n",
    "            \"p_gt0\": p_gt0,\n",
    "            \"hdi95_low\": float(hdi95[0]), \"hdi95_high\": float(hdi95[1]),\n",
    "        })\n",
    "        print(f\"[{label}] {metric}: cached -> {fname}\")\n",
    "\n",
    "# write a small index for convenience\n",
    "pd.DataFrame(index_rows).to_csv(f\"{CACHE_DIR}/index.csv\", index=False)\n",
    "print(\"Done. Cache written to\", CACHE_DIR)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14a95e99",
   "metadata": {},
   "source": [
    "# Create table of posterior distributions for each metric and case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "4edf2ade",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(a)] wrote bayes_post_a.tex\n",
      "[(b)] wrote bayes_post_b.tex\n",
      "[(c1)] wrote bayes_post_c1.tex\n",
      "[(c2)] wrote bayes_post_c2.tex\n"
     ]
    }
   ],
   "source": [
    "CACHE_DIR = \"bayes_cache\"\n",
    "metrics_order = [\"Accuracy\",\"Precision\",\"Sensitivity\",\"Specificity\",\"F1 Score\",\"MCC\"]\n",
    "\n",
    "def load_one(npz_path):\n",
    "    d = np.load(npz_path, allow_pickle=True)\n",
    "    rec = {k: (d[k].item() if d[k].shape == () else d[k]) for k in d.files}\n",
    "    # fallbacks if old cache (without diagnostics)\n",
    "    for miss in [\"hdi90_low\",\"hdi90_high\",\"rhat_b0\",\"rhat_b1\",\"ess_bulk_b0\",\"ess_bulk_b1\"]:\n",
    "        if miss not in rec:\n",
    "            # compute 90% HDI from draws if needed\n",
    "            if miss in [\"hdi90_low\",\"hdi90_high\"]:\n",
    "                lo, hi = az.hdi(rec[\"beta1\"], hdi_prob=0.90)\n",
    "                rec[\"hdi90_low\"], rec[\"hdi90_high\"] = float(lo), float(hi)\n",
    "            else:\n",
    "                rec[miss] = np.nan\n",
    "    return rec\n",
    "\n",
    "def fmt_mean_sd(arr):\n",
    "    return f\"{np.mean(arr):.4f} ({np.std(arr, ddof=1):.4f})\"\n",
    "\n",
    "def fmt_prob(p):\n",
    "    return f\"{p:.2f}\"\n",
    "\n",
    "def build_tables():\n",
    "    files = glob.glob(os.path.join(CACHE_DIR, \"post_*.npz\"))\n",
    "    if not files:\n",
    "        raise SystemExit(\"No cached files found in bayes_cache/\")\n",
    "    rows = []\n",
    "    for f in files:\n",
    "        rec = load_one(f)\n",
    "        lab = str(rec[\"label\"])\n",
    "        metric = str(rec[\"metric\"])\n",
    "        b0 = rec[\"beta0\"]; b1 = rec[\"beta1\"]\n",
    "        rows.append({\n",
    "            \"label\": lab,\n",
    "            \"metric\": metric,\n",
    "            \"beta0_mean_sd\": fmt_mean_sd(b0),\n",
    "            \"beta1_mean_sd\": fmt_mean_sd(b1),\n",
    "            \"p_beta1_gt0\": rec[\"p_gt0\"],\n",
    "            \"hdi90_low\": rec[\"hdi90_low\"],\n",
    "            \"hdi90_high\": rec[\"hdi90_high\"],\n",
    "            \"ess_bulk_b0\": rec[\"ess_bulk_b0\"],\n",
    "            \"ess_bulk_b1\": rec[\"ess_bulk_b1\"],\n",
    "            \"rhat_b0\": rec[\"rhat_b0\"],\n",
    "            \"rhat_b1\": rec[\"rhat_b1\"],\n",
    "            \"model_1\": rec[\"model_1\"],\n",
    "            \"model_2\": rec[\"model_2\"],\n",
    "        })\n",
    "    df = pd.DataFrame(rows)\n",
    "\n",
    "    # One LaTeX table per comparison label\n",
    "    for lab, g in df.groupby(\"label\", sort=False):\n",
    "        # keep metric order\n",
    "        g = g.set_index(\"metric\").reindex(metrics_order).reset_index().dropna(subset=[\"beta0_mean_sd\"])\n",
    "\n",
    "        #cap = (f\"Bayesian GLM posterior summaries for comparison {lab}: \"\n",
    "        #       f\"Model 1 = {g['model_1'].iloc[0]}, Model 2 = {g['model_2'].iloc[0]}. \"\n",
    "        #       \"Means and SDs are over posterior draws; HDI is 90%; ESS is bulk ESS; \"\n",
    "        #       \"R̂ is split-\\u0052hat.\")\n",
    "        \n",
    "        cap = (f\"Bayesian GLM posterior summaries for comparison {lab}: \"\n",
    "               f\"Model 1 = {g['model_1'].iloc[0]}, Model 2 = {g['model_2'].iloc[0]}.\")\n",
    "        label_tex = f\"tab:bayes_post_{lab.strip('()')}\"\n",
    "        colspec = \"l l l r r r r\"\n",
    "\n",
    "        # assemble a compact table with β0 and β1 in separate columns\n",
    "        out = []\n",
    "        out.append(\"\\\\begin{table}[t]\")\n",
    "        out.append(\"\\\\footnotesize\")\n",
    "        out.append(f\"\\\\caption{{{cap}}}\")\n",
    "        out.append(f\"\\\\label{{{label_tex}}}\")\n",
    "        out.append(f\"\\\\begin{{tabular}}{{l l l r r r r}}\")\n",
    "        out.append(\"\\\\toprule\")\n",
    "        out.append(\"Metric & $\\\\beta_0$ mean (SD) & $\\\\beta_1$ mean (SD) & $P(\\\\beta_1>0)$ & ESS$_{\\\\beta_0}$ & ESS$_{\\\\beta_1}$ & $\\\\hat R_{\\\\beta_0}$ / $\\\\hat R_{\\\\beta_1}$\\\\\\\\\")\n",
    "        out.append(\"\\\\midrule\")\n",
    "\n",
    "        for _, r in g.iterrows():\n",
    "            #hdi = f\"{r['hdi90_low']:.4f}, {r['hdi90_high']:.4f}\"\n",
    "            ess_b0 = \"—\" if np.isnan(r[\"ess_bulk_b0\"]) else f\"{r['ess_bulk_b0']:.0f}\"\n",
    "            ess_b1 = \"—\" if np.isnan(r[\"ess_bulk_b1\"]) else f\"{r['ess_bulk_b1']:.0f}\"\n",
    "            rhat_b0 = \"—\" if np.isnan(r[\"rhat_b0\"]) else f\"{r['rhat_b0']:.3f}\"\n",
    "            rhat_b1 = \"—\" if np.isnan(r[\"rhat_b1\"]) else f\"{r['rhat_b1']:.3f}\"\n",
    "            out.append(\n",
    "                f\"{r['metric']} & {r['beta0_mean_sd']} & {r['beta1_mean_sd']} & {fmt_prob(r['p_beta1_gt0'])} & \"\n",
    "                f\"{ess_b0} & {ess_b1} & {rhat_b0} / {rhat_b1} \\\\\\\\\"\n",
    "            )\n",
    "\n",
    "        out.append(\"\\\\bottomrule\")\n",
    "        out.append(\"\\\\end{tabular}\")\n",
    "        out.append(\"\\\\end{table}\")\n",
    "\n",
    "        tex_path = f\"bayes_post_{lab.strip('()')}.tex\"\n",
    "        with open(tex_path, \"w\", encoding=\"utf-8\") as fh:\n",
    "            fh.write(\"\\n\".join(out))\n",
    "        print(f\"[{lab}] wrote {tex_path}\")\n",
    "\n",
    "build_tables()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "15fb57ab",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(a)] wrote bayes_violin_beta0_a.pdf\n",
      "[(a)] wrote bayes_violin_beta1_a.pdf\n",
      "[(b)] wrote bayes_violin_beta0_b.pdf\n",
      "[(b)] wrote bayes_violin_beta1_b.pdf\n",
      "[(c1)] wrote bayes_violin_beta0_c1.pdf\n",
      "[(c1)] wrote bayes_violin_beta1_c1.pdf\n",
      "[(c2)] wrote bayes_violin_beta0_c2.pdf\n",
      "[(c2)] wrote bayes_violin_beta1_c2.pdf\n"
     ]
    }
   ],
   "source": [
    "CACHE_DIR = \"bayes_cache\"\n",
    "PR_BLUE = 0.95   # blue if P(β1>0) >= 0.95  -> model 2 better\n",
    "PR_RED  = 0.05   # red  if P(β1>0) <= 0.05  -> model 2 worse\n",
    "\n",
    "metrics = [\"Accuracy\", \"Precision\", \"Sensitivity\", \"Specificity\", \"F1 Score\", \"MCC\"]\n",
    "comparisons = [\n",
    "    (\"(a)\", \"basic_no_fft\", \"cdd_fft\"),\n",
    "    (\"(b)\", \"basic_no_fft\", \"cdd_no_fft\"),\n",
    "    (\"(c1)\", \"basic_no_fft\", \"basic_fft\"),\n",
    "    (\"(c2)\", \"cdd_no_fft\", \"cdd_fft\"),\n",
    "]\n",
    "\n",
    "# Colors\n",
    "GREY = \"#808080\"\n",
    "BLUE = \"#1f77b4\"\n",
    "RED  = \"#d62728\"\n",
    "\n",
    "def color_for_p(p_gt0: float) -> str:\n",
    "    if p_gt0 >= PR_BLUE:\n",
    "        return BLUE\n",
    "    elif p_gt0 <= PR_RED:\n",
    "        return RED\n",
    "    return GREY\n",
    "\n",
    "def load_cached(label: str, metric: str):\n",
    "    \"\"\"Load one cached npz produced by your sampling step.\"\"\"\n",
    "    lab = label.strip(\"()\")\n",
    "    f = os.path.join(CACHE_DIR, f\"post_{lab}__{metric.replace(' ', '_')}.npz\")\n",
    "    if not os.path.exists(f):\n",
    "        return None\n",
    "    data = np.load(f, allow_pickle=True)\n",
    "    return {\n",
    "        \"beta0\": data[\"beta0\"],\n",
    "        \"beta1\": data[\"beta1\"],\n",
    "        \"p_gt0\": float(data[\"p_gt0\"]),\n",
    "        \"hdi95_low\": float(data[\"hdi95_low\"]),\n",
    "        \"hdi95_high\": float(data[\"hdi95_high\"]),\n",
    "        \"label\": str(data[\"label\"]),\n",
    "        \"model_1\": str(data[\"model_1\"]),\n",
    "        \"model_2\": str(data[\"model_2\"]),\n",
    "        \"metric\": str(data[\"metric\"]),\n",
    "    }\n",
    "\n",
    "def _add_mean_bar(ax, x, samples, width=0.4, **kwargs):\n",
    "    \"\"\"\n",
    "    Draw a short horizontal bar at the mean of `samples`\n",
    "    centred on x, with total length `width`.\n",
    "    \"\"\"\n",
    "    m = float(np.mean(samples))\n",
    "    ax.hlines(m, x - width/2, x + width/2, **{\"colors\": \"black\", \"linewidth\": 2, **kwargs})\n",
    "\n",
    "def _style_grid(ax):\n",
    "    ax.grid(True, which=\"major\", linestyle=\":\", linewidth=0.7, alpha=0.6)\n",
    "\n",
    "def plot_beta0(entries, label, model_1, model_2):\n",
    "    \"\"\"One figure: β0 violins (grey) for all metrics in this comparison.\"\"\"\n",
    "    nM = len(entries)\n",
    "    fig, ax = plt.subplots(figsize=(1.0 * nM, 4.0), constrained_layout=True)\n",
    "\n",
    "    x_positions = np.arange(nM)\n",
    "    for i, r in enumerate(entries):\n",
    "        parts = ax.violinplot(r[\"beta0\"], positions=[i], widths=0.8, showextrema=False)\n",
    "        for pc in parts[\"bodies\"]:\n",
    "            pc.set_facecolor(\"#bbbbbb\")\n",
    "            pc.set_edgecolor(\"black\")\n",
    "            pc.set_alpha(1.0)\n",
    "        # Mean bar inside each violin\n",
    "        _add_mean_bar(ax, i, r[\"beta0\"])\n",
    "\n",
    "    # Axis formatting\n",
    "    ax.set_xticks(x_positions)\n",
    "    ax.set_xticklabels([r[\"metric\"] for r in entries], rotation=45, ha=\"right\")\n",
    "    ax.set_ylabel(\"Posterior draws (β₀)\")\n",
    "\n",
    "    # (e) Fixed y-axis and ticks\n",
    "    ax.set_ylim(0.55, 0.99)\n",
    "    ax.set_yticks(np.arange(0.55, 0.99 + 1e-9, 0.1))\n",
    "\n",
    "    # Grid\n",
    "    _style_grid(ax)\n",
    "\n",
    "    #title = f\"Bayesian GLM posteriors — β₀ only — case {label} — Model 1={model_1}, Model 2={model_2}\"\n",
    "    title = f\"Posteriors for β₀ — case {label} — Model 1={model_1}, Model 2={model_2}\"\n",
    "    #subtitle = f\"Baseline mean per metric; Model 1={model_1}, Model 2={model_2}\"\n",
    "    #ax.set_title(title + \"\\n\" + subtitle, fontsize=10)\n",
    "    ax.set_title(title, fontsize=10)\n",
    "\n",
    "    # Legend\n",
    "    ax.legend(\n",
    "        handles=[mpatches.Patch(facecolor=\"#bbbbbb\", edgecolor=\"black\", label=\"β₀ (baseline mean)\")\n",
    "        ],\n",
    "        loc=\"lower left\", fontsize=10, frameon=False\n",
    "    )\n",
    "    out = f\"bayes_violin_beta0_{label.strip('()')}.pdf\"\n",
    "    fig.savefig(out, dpi=300, bbox_inches=\"tight\")\n",
    "    plt.close(fig)\n",
    "    print(f\"[{label}] wrote {out}\")\n",
    "\n",
    "def plot_beta1(entries, label, model_1, model_2):\n",
    "    \"\"\"One figure: β1 violins (colored by P(β1>0)) for all metrics in this comparison.\"\"\"\n",
    "    nM = len(entries)\n",
    "    fig, ax = plt.subplots(figsize=(1.0 * nM, 4.2), constrained_layout=True)\n",
    "\n",
    "    x_positions = np.arange(nM)\n",
    "    for i, r in enumerate(entries):\n",
    "        col = color_for_p(r[\"p_gt0\"])\n",
    "        parts = ax.violinplot(r[\"beta1\"], positions=[i], widths=0.8, showextrema=False)\n",
    "        for pc in parts[\"bodies\"]:\n",
    "            pc.set_facecolor(col)\n",
    "            pc.set_edgecolor(\"black\")\n",
    "            pc.set_alpha(1.0)\n",
    "        # Mean bar inside each violin\n",
    "        _add_mean_bar(ax, i, r[\"beta1\"])\n",
    "\n",
    "    # zero reference for β1\n",
    "    ax.axhline(0.0, color=\"black\", linewidth=0.8, alpha=0.7)\n",
    "\n",
    "    # Axis formatting\n",
    "    ax.set_xticks(x_positions)\n",
    "    ax.set_xticklabels([r[\"metric\"] for r in entries], rotation=45, ha=\"right\")\n",
    "    ax.set_ylabel(\"Posterior draws (β₁)\")\n",
    "\n",
    "    # (d) Fixed y-axis and ticks\n",
    "    ax.set_ylim(-0.065, 0.135)\n",
    "    ax.set_yticks(np.arange(-0.06, 0.13 + 1e-9, 0.02))\n",
    "\n",
    "    # Grid\n",
    "    _style_grid(ax)\n",
    "\n",
    "    #title = f\"Bayesian GLM posteriors — β₁ only — case {label} — Model 1={model_1}, Model 2={model_2}\"\n",
    "    title = f\"Posteriors for β₁ — case {label} — Model 1={model_1}, Model 2={model_2}\"\n",
    "    #subtitle = (\n",
    "    #    f\"Blue if P(β₁>0)≥{PR_BLUE:.2f} (Model 2 better), \"\n",
    "    #    f\"red if ≤{PR_RED:.2f} (Model 2 worse), grey otherwise; \"\n",
    "    #    f\"Model 1={model_1}, Model 2={model_2}\"\n",
    "    #)\n",
    "    #ax.set_title(title + \"\\n\" + subtitle, fontsize=10)\n",
    "    ax.set_title(title, fontsize=10)\n",
    "\n",
    "    # (a) Removed the per-column P(β₁>0)=... headings\n",
    "\n",
    "    # Legend\n",
    "    leg_items = [\n",
    "        mpatches.Patch(facecolor=BLUE, edgecolor=\"black\", label=f\"Model 2 better (P(β₁>0)≥{PR_BLUE:.2f})\"),\n",
    "        mpatches.Patch(facecolor=RED,  edgecolor=\"black\", label=f\"Model 2 worse (P(β₁>0)≤{PR_RED:.2f})\"),\n",
    "        mpatches.Patch(facecolor=GREY, edgecolor=\"black\", label=\"no significant difference\"),\n",
    "    ]\n",
    "    ax.legend(handles=leg_items, loc=\"upper left\", fontsize=10, frameon=False)\n",
    "\n",
    "    out = f\"bayes_violin_beta1_{label.strip('()')}.pdf\"\n",
    "    fig.savefig(out, dpi=300, bbox_inches=\"tight\")\n",
    "    plt.close(fig)\n",
    "    print(f\"[{label}] wrote {out}\")\n",
    "\n",
    "# ---------------- Main ----------------\n",
    "if __name__ == \"__main__\":\n",
    "    for label, model_1, model_2 in comparisons:\n",
    "        # Load all metrics available for this comparison\n",
    "        entries = []\n",
    "        for metric in metrics:\n",
    "            rec = load_cached(label, metric)\n",
    "            if rec is not None:\n",
    "                entries.append(rec)\n",
    "\n",
    "        if not entries:\n",
    "            print(f\"[{label}] no cached files found in {CACHE_DIR} — skip\")\n",
    "            continue\n",
    "\n",
    "        # keep original metric order\n",
    "        entries = sorted(entries, key=lambda r: metrics.index(r[\"metric\"]))\n",
    "\n",
    "        plot_beta0(entries, label, model_1, model_2)\n",
    "        plot_beta1(entries, label, model_1, model_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f6a8b6c",
   "metadata": {},
   "source": [
    "# Hypothesis testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "4676b43c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# -------- Load and prepare data --------\n",
    "df = pd.read_csv(\"results_metrics.csv\")\n",
    "\n",
    "# Filter to External dataset\n",
    "df = df[df[\"Dataset\"] == \"External\"].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "cbc0090c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Corrected mapping for model names\n",
    "df[\"Model\"] = df[\"Model\"].replace({\n",
    "    \"Basic Model No FFT\": \"basic_no_fft\",\n",
    "    \"Basic Model With FFT\": \"basic_fft\",\n",
    "    \"CDD-based Model No FFT\": \"cdd_no_fft\",\n",
    "    \"CDD-based Model With FFT\": \"cdd_fft\",\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "76aa8e9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Metrics and comparisons\n",
    "metrics = [\"Accuracy\", \"Precision\", \"Sensitivity\", \"Specificity\", \"F1 Score\", \"MCC\"]\n",
    "comparisons = [\n",
    "    (\"(a)\", \"cdd_fft\", \"basic_no_fft\"),\n",
    "    (\"(b)\", \"cdd_no_fft\", \"basic_no_fft\"),\n",
    "    (\"(c1)\", \"basic_fft\", \"basic_no_fft\"),\n",
    "    (\"(c2)\", \"cdd_fft\", \"cdd_no_fft\"),\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "1f1b4bb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---------- Helpers ----------\n",
    "def hedges_g(x, y):\n",
    "    \"\"\"Hedges' g (independent samples).\"\"\"\n",
    "    n1, n2 = len(x), len(y)\n",
    "    if n1 < 2 or n2 < 2:\n",
    "        return np.nan\n",
    "    s1, s2 = np.var(x, ddof=1), np.var(y, ddof=1)\n",
    "    dfp = n1 + n2 - 2\n",
    "    if dfp <= 0:\n",
    "        return np.nan\n",
    "    sp2 = ((n1 - 1)*s1 + (n2 - 1)*s2) / dfp\n",
    "    if sp2 <= 0:\n",
    "        return np.nan\n",
    "    sp = np.sqrt(sp2)\n",
    "    d = (np.mean(x) - np.mean(y)) / sp\n",
    "    J = 1 - (3 / (4*(n1 + n2) - 9)) if (n1 + n2) > 2 else 1.0\n",
    "    return J * d\n",
    "\n",
    "def mean_sd_str(a):\n",
    "    return f\"{np.mean(a):.4f} ({np.std(a, ddof=1):.4f})\"\n",
    "\n",
    "def latex_sci(p, sig=2, min_exp=-308):\n",
    "    \"\"\"Format p as x·10^{y} in LaTeX math mode. Handle underflow gracefully.\"\"\"\n",
    "    if p <= 0 or not np.isfinite(p):\n",
    "        return rf\"$< 1\\cdot10^{{{min_exp}}}$\"\n",
    "    exp = int(np.floor(np.log10(p)))\n",
    "    mant = p / (10**exp)\n",
    "    return rf\"${mant:.{sig}f}\\cdot10^{{{exp}}}$\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "65e03b41",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(a)] wrote welch_table_a.tex\n",
      "[(b)] wrote welch_table_b.tex\n",
      "[(c1)] wrote welch_table_c1.tex\n",
      "[(c2)] wrote welch_table_c2.tex\n"
     ]
    }
   ],
   "source": [
    "# ---------- Make one table per comparison ----------\n",
    "def write_table(label, model1, model2):\n",
    "    rows = []\n",
    "    for metric in metrics:\n",
    "        x = df.loc[df[\"Model\"] == model1, metric].dropna().values\n",
    "        y = df.loc[df[\"Model\"] == model2, metric].dropna().values\n",
    "        if len(x) < 2 or len(y) < 2:\n",
    "            continue\n",
    "        # Welch’s t-test (two-sided)\n",
    "        t_stat, p_val = stats.ttest_ind(x, y, equal_var=False)\n",
    "\n",
    "        rows.append({\n",
    "            \"Metric\": metric,\n",
    "            \"Model 1 mean (SD)\": mean_sd_str(x),\n",
    "            \"Model 2 mean (SD)\": mean_sd_str(y),\n",
    "            \"Difference\": np.mean(x) - np.mean(y),      # M1 - M2\n",
    "            \"t\": t_stat,\n",
    "            \"Hedges g\": hedges_g(x, y),\n",
    "            \"p\": latex_sci(float(p_val), sig=2)         # x·10^{y}\n",
    "        })\n",
    "\n",
    "    if not rows:\n",
    "        print(f\"[{label}] No valid metrics for {model1} vs {model2}\")\n",
    "        return\n",
    "\n",
    "    res = pd.DataFrame(rows).sort_values(\"Metric\").reset_index(drop=True)\n",
    "\n",
    "    # Build LaTeX manually to inject \\footnotesize and booktabs lines\n",
    "    colspec = \"l l l r r r r\"\n",
    "    header = [\"Metric\", \"Model 1 mean (SD)\", \"Model 2 mean (SD)\", \"Difference\", \"t\", \"Hedges g\", \"p\"]\n",
    "    caption = (\n",
    "        f\"Welch's t-test on the External dataset for comparison {label}: \"\n",
    "        f\"Model 1 = {model1.replace('_','\\\\_')} vs. Model 2 = {model2.replace('_','\\\\_')}, \"\n",
    "        \"SD: Standard Deviation.\"\n",
    "    )\n",
    "    label_tex = f\"tab:welch_external_{label.strip('()')}\"\n",
    "\n",
    "    lines = []\n",
    "    lines.append(\"\\\\begin{table}\")\n",
    "    lines.append(\"\\\\footnotesize\")\n",
    "    lines.append(f\"\\\\caption{{{caption}}}\")\n",
    "    lines.append(f\"\\\\label{{{label_tex}}}\")\n",
    "    lines.append(f\"\\\\begin{{tabular}}{{{colspec}}}\")\n",
    "    lines.append(\"\\\\toprule\")\n",
    "    lines.append(\" & \".join(header) + \" \\\\\\\\\")\n",
    "    lines.append(\"\\\\midrule\")\n",
    "\n",
    "    for _, r in res.iterrows():\n",
    "        line = \" & \".join([\n",
    "            r[\"Metric\"],\n",
    "            r[\"Model 1 mean (SD)\"],\n",
    "            r[\"Model 2 mean (SD)\"],\n",
    "            f\"{r['Difference']:.4f}\",\n",
    "            f\"{r['t']:.4f}\",\n",
    "            f\"{r['Hedges g']:.4f}\" if np.isfinite(r[\"Hedges g\"]) else \"\",\n",
    "            r[\"p\"],\n",
    "        ]) + \" \\\\\\\\\"\n",
    "        lines.append(line)\n",
    "\n",
    "    lines.append(\"\\\\bottomrule\")\n",
    "    lines.append(\"\\\\end{tabular}\")\n",
    "    lines.append(\"\\\\end{table}\")\n",
    "\n",
    "    out_path = f\"welch_table_{label.strip('()')}.tex\"\n",
    "    with open(out_path, \"w\", encoding=\"utf-8\") as f:\n",
    "        f.write(\"\\n\".join(lines))\n",
    "    print(f\"[{label}] wrote {out_path}\")\n",
    "\n",
    "for label, m1, m2 in comparisons:\n",
    "    write_table(label, m1, m2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0ee1f79",
   "metadata": {},
   "source": [
    "# Old Bayesian analysis code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "a0bb5d3b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [mu1, mu2, sigma]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2ece9e53e0bc4ddc94b9862e3b760992",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 21 seconds.\n"
     ]
    }
   ],
   "source": [
    "# Extract the metric values for the two models\n",
    "y1 = df[df[\"Model\"] == model_1][metric].values\n",
    "y2 = df[df[\"Model\"] == model_2][metric].values\n",
    "\n",
    "# Aaaaaaaaand run the Bayesian model\n",
    "with pm.Model() as model:\n",
    "    mu1 = pm.Normal(\"mu1\", mu=0.5, sigma=1)\n",
    "    mu2 = pm.Normal(\"mu2\", mu=0.5, sigma=1)\n",
    "    sigma = pm.Exponential(\"sigma\", lam=10)\n",
    "    #sigma = pm.HalfNormal(\"sigma\", sigma=0.1)\n",
    "\n",
    "    y1_obs = pm.Normal(\"y1_obs\", mu=mu1, sigma=sigma, observed=y1)\n",
    "    y2_obs = pm.Normal(\"y2_obs\", mu=mu2, sigma=sigma, observed=y2)\n",
    "\n",
    "    delta = pm.Deterministic(\"delta\", mu2 - mu1)\n",
    "\n",
    "    trace = pm.sample(2000, tune=1000, target_accept=0.95, return_inferencedata=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "c0e90cac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.microsoft.datawrangler.viewer.v0+json": {
       "columns": [
        {
         "name": "index",
         "rawType": "object",
         "type": "string"
        },
        {
         "name": "mean",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "sd",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "hdi_5%",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "hdi_95%",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "mcse_mean",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "mcse_sd",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "ess_bulk",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "ess_tail",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "r_hat",
         "rawType": "float64",
         "type": "float"
        }
       ],
       "ref": "5700e636-a8cd-4193-bb7f-b35623f7e9f1",
       "rows": [
        [
         "mu1",
         "0.631",
         "0.003",
         "0.626",
         "0.635",
         "0.0",
         "0.0",
         "6523.0",
         "5009.0",
         "1.0"
        ],
        [
         "mu2",
         "0.73",
         "0.003",
         "0.725",
         "0.734",
         "0.0",
         "0.0",
         "6850.0",
         "4607.0",
         "1.0"
        ],
        [
         "delta",
         "0.099",
         "0.004",
         "0.093",
         "0.105",
         "0.0",
         "0.0",
         "6690.0",
         "5092.0",
         "1.0"
        ]
       ],
       "shape": {
        "columns": 9,
        "rows": 3
       }
      },
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_5%</th>\n",
       "      <th>hdi_95%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>mu1</th>\n",
       "      <td>0.631</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.626</td>\n",
       "      <td>0.635</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>6523.0</td>\n",
       "      <td>5009.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mu2</th>\n",
       "      <td>0.730</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.725</td>\n",
       "      <td>0.734</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>6850.0</td>\n",
       "      <td>4607.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>delta</th>\n",
       "      <td>0.099</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.093</td>\n",
       "      <td>0.105</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>6690.0</td>\n",
       "      <td>5092.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        mean     sd  hdi_5%  hdi_95%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\n",
       "mu1    0.631  0.003   0.626    0.635        0.0      0.0    6523.0    5009.0   \n",
       "mu2    0.730  0.003   0.725    0.734        0.0      0.0    6850.0    4607.0   \n",
       "delta  0.099  0.004   0.093    0.105        0.0      0.0    6690.0    5092.0   \n",
       "\n",
       "       r_hat  \n",
       "mu1      1.0  \n",
       "mu2      1.0  \n",
       "delta    1.0  "
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Print summary and plot posterior distributions\n",
    "az.summary(trace, var_names=[\"mu1\", \"mu2\", \"delta\"], hdi_prob=0.90)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "4487160e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAG3CAYAAAAU+jfPAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZzJJREFUeJzt3Qd4U9X7B/C3e+8FpZuyR9l7imwFURAVZboRwb1+TtwKDlTc4B+3KCrIUhFkz7Ipo7RQRumke/f+n/ekNyZp2qZtkpvx/TxPSbi5uTm5ueO957znXAdJkiQCAAAAMBNHc30QAAAAAEPwAQAAAGaF4AMAAADMCsEHAAAAmBWCDwAAADArBB8AAABgVgg+AAAAwKwQfAAAAIBZIfgAAAAAs0LwYUWWL19ODg4OlJqaSpZk7969NGDAAPLy8hLlO3jwINmqYcOGiT8Z/xb8nfm30bR+/Xrq1q0bubu7i9evXr0qpq9YsYLat29PLi4u5O/vb/by27qYmBiaOXOmyfalurb1un5vU5C3ubfffpvMidcrr19bgX1UWTYdfMgHGPmPN7K2bdvSAw88QFeuXDH65xUXF9MLL7xAmzdvJntRUVFBU6ZMoZycHHrnnXfEjhsdHa13Xl4v8m/x9ddf651n4MCB4vXOnTvXeq2qqoqWLVsmTv6BgYHk5uYmDoazZs2iffv21Zo/OTmZ7rnnHoqLixO/va+vr1j+e++9RyUlJWQq2dnZdPPNN5OHhwd9+OGHYp3wySopKUkcwFu3bk2fffYZffrppyYrA5hvW6/r93711Vfp119/VbrY0Mx91B6P62Yh2bBly5bxfWukl156SVqxYoX02WefSTNmzJAcHR2l2NhYqaioyKifl5mZKT7v+eefl0yhsrJSKikpkaqrqyVLceLECfGded025J9//hHzuru7S2PHjq31ekpKivr1Tp06ab1WXFwsjRkzRrw+ZMgQ6a233pK++OIL6dlnn5XatWsnOTg4SGlpaer516xZI3l4eEj+/v7Sgw8+KH366afSBx98IN1yyy2Si4uLdNdddzXp+w4dOlT8yfi34N+EfxvZunXrRDn//PNPrfcuXbpUTD99+nSTPhsaFh0dLfZxQ48NvM01d1uv6/f28vIyqCyNJe8nvA+YU3l5uVRaWirZgsbso6Y+rtsrZ7IDY8eOpV69eonnd955JwUFBdHixYvpt99+o1tvvZUsXVFRkYjKnZycxJ+xcETv6enZrGVkZGSIx8ZUT44bN45+//13ysrKouDgYPX0b7/9lsLCwqhNmzaUm5ur9Z7HHntMVJPyFeeCBQu0Xnv++efFdFlKSgrdcsst4qp006ZN1LJlS/Vrc+fOpTNnztAff/xBxiDXqBmyTpqyrgzdNsD0zPm7WiJuhrAV9v5bWgTJhslXN3v37tWazlfFPP2VV14R/6+oqBC1I3FxcZKrq6u4enrqqadqRfm8nFGjRklBQUHi6jwmJkaaNWuW1tWI7p9mtMxXTjfddJMUEBAgubm5ST179pR+++03vWXevHmzdN9990khISHi6r2+q7UPP/xQ6tixoyh7y5Ytpfvvv1/Kzc3Vmoev1rk2Yd++fdLgwYNFrcD8+fPrXX9///23NGjQIMnT01Py8/OTJkyYIB0/flz9Ol/V6X5fzVqBumo+vvrqK3FV+NFHH2m9zuWbN2+euqwyrtFwdnaWRo4cKRni3nvvFZ+zfft2qTk++eQTsU3wb927d2/p33//rVXzIf/u/Nswfk13nfB64m2qvm1j7dq16nXt7e0tjRs3Tjp69KhWeXg5vN7OnDkjao54vokTJ4rXqqqqpHfeeUdsB7xthYaGSnfffbeUk5OjtQwux/jx46WtW7eK78Tzci0g/ya6eBtasGCBeA9vW61atZLuuOMOcSUo433kueeek1q3bi3miYiIkB577LFa+w6/h7d/Q2sbuaaSyyfXXvE2u2HDBq0ap4ULF4oy8TzDhg0T60tfzQdPHz58uPgdeX5+H9eaNabmo65tva7fW9+xwFi1IJo1H4sXL5aioqLEd+MawSNHjmjNe+jQIfG5/Bvzbx0WFiaOWVlZWVrz5efni+OB/Fvzcefaa6+V9u/fr7UO+HVNvN29++67UufOncXyg4ODpdGjR9c65tZH3q4vXLggtmd+zst55JFHtGoUWWFhofTwww+L7YzL2bZtW7EeGlMb3Jh9tK7fErUgzWcXNR/6cgEY14DItSFfffUVTZ48mR555BHavXs3vfbaa3TixAlatWqVOiIeNWoUhYSE0JNPPikiY078+uWXX8TrPH3p0qV033330aRJk+jGG28U07t27Soejx07JvINWrVqJd7PV6s//vgj3XDDDfTzzz+L92i6//77xTKfe+45cXVbF26LfPHFF+naa68Vn33y5ElRDk6M2759u9bVCrdzci0Q1wrcfvvtopahLn/99ZeYl/Ml+DM4R2LJkiXiOxw4cEDkWnA+BX8fbtt+8MEHqXfv3vUuU8a1LRMnTqTvvvtOlJkdOnRIrKPPP/+cDh8+rDX/unXrqLKyku644w4yxOrVq0W5OTGwqb744gvx/XgZXNNy9uxZmjBhgsg1iYyMrPN9zzzzDLVr1060Fb/00ksUGxsr2o/5d/6///s/sT3x7+Pt7a3eNri9ecaMGTR69Gh64403RI0UzzNo0CBKTEzUSvLj9cDz8WuccCjXXHFZOceJ81/4t+Danw8++EC8X3c74Jof3tbnzJkjPvfLL78U7dw9e/akTp06iXkKCwtp8ODBYh+YPXs29ejRQ9RUcY3VhQsXRI1VdXW1WCfbtm2ju+++mzp06EBHjhwRtVCnTp3SynfgsvB2+s8//2gl7OrD8/E2x+ue16Grq6vYJ7kWi/dBxvvFyy+/LGrR+I+3SX6tvLxca1np6ek0fPhwsd7k/Y5/G27rb4z6tnV9vzfvj3xc6dOnj1g3jKcbE29PBQUFojavtLRU5DJdc8014jeQy/bnn3+KbZe3ixYtWoh9jMvKj7t27RI1d+zee++llStXiny4jh07imMF/678+/NvXxfehni742MFf19ez1u3bhXLlmubDcH5XLxd9+3bV2zXfPxZtGiRWGfyMUKSJLG98TbEn8vJohs2bBC1ohcvXtSq/axPY/bRLl26UL9+/eo8rkMzSDZMrin466+/xJUXX0F///33ouaCr5Y40j548KCY584779R676OPPiqmb9q0Sfx/1apVemtRNNXXNjhixAipS5cuWleEHK0PGDBAatOmTa0y81WwbtSvW/ORkZEhon+ujeErEBnnNvB8X375Za1o/+OPPzZo3XXr1k1cPWdnZ2tdRXG+zPTp02vVZvz0008NLlNzXq594jyN8+fPi9f4aplrGeSyatZ8PPTQQ+J9iYmJDX5GXl6emFeuEWhq2zZ/d14HZWVl6umcN6Jbu6Nb81FfjRtvFzxds+agoKBAXNnr5qCkp6eL2ibN6fJV2JNPPqk1L9di8PRvvvlGa/r69etrTZev7rgWR8bbEV+18pWmjGszeL5ffvml1vqRrzK5doK3B/58TbyN6dY8yd+dt4H6cFs7L3PSpEla27Tm58rbPdfgaF7xPv3007VqGLjmhqft3r1b6/vyum1szkdd23pdv7epcz7kY5iMvyNP5/1FM1dK13fffVdrG+D1MXfu3Ho/V7fmg4+NvBzOqdLVmJoIebvm2mdN3bt3F7XDsl9//VXM9/LLL2vNN3nyZHEs4RpBQzVmH0XOh2nYdG8XGV+FcC0CX7HyVT9HtBzd8pXM2rVrxTwPP/yw1nu4BoTJuQFyG+CaNWtE1ntjcHY8X7VxdjVfqfAVJP/x1QVH+6dPnxaRu6a77rqrwfwOvjrgKz2+Mnd0dNR6L/fs0M1r4N4hfAXUkMuXL4suhHw1zFf6Mo72R44cqV5nzcFXqbzs77//XlzR8GNd+Tf5+fni0cfHp8HlNmbeunDPGa7p4qtBvuqW8frw8/MjY+IrU+7ix99d3i74j397vgrkqzxd8pWg7KeffhLl4t9Gcxlck8Hbuu4y+MqWazVkvG/wlSBfIcu4Ni4hIaFWjRyTr5b5c7m2g7slan4uX30zzc/lmgz+nRuq9eDaEq5R4ZoNzW1a83Pl7X7evHnqaUw3F4jxtspXrlwDofl9p02bRtaOr9T5GCbj78jbjOb+qVnDw7Uj/Pvw+mBcWyTj4xvXLl26dMngz+dthNc/51zp0vxdDMX7mybeRjW3Sf5evF9wzZPusZq3La4hBethF80u3JWKu9g6OzuL6kg+0MoHtnPnzonn8fHxWu/hKkreIfl1NnToULrppptElTBX7/FBlHf+2267TZzU68PV3LxzPPvss+JPHz7ZaR5IuCqwIXLZ+Pto4hMmNzvIr8t4+Zon08Yul/HJhqs6m5voyM0A3G2Rk0z5oJmWlibWpT4cSDEO3BrSmHkb+v6c+KpbZl6vxsSBJ5NP2HV9HxlvwxEREbWWkZeXR6GhoXqXISfRyaKiomrNExAQoJXky02TvL03VHaulueTuSGfawj+XN4fOUBq7O/D5eDvoTsvn5B16du2lcTNQ5o4mGyoaUj3+zM+znFzruaFDx+zOLjX/T14m5G9+eabogmOL9A4aOWmrOnTp9e7vfNvFR4ernWB0lSctK27Heluk/xb8ufpXljwMUl+HayHXQQffHJrqP2xoUidX+c2UW7L5JwCPgFzWzi3S/I0vsKsC1/JsUcffVTUdOijG/w0tk3aEKZYZnNwsPHxxx+Lq2K+yq7rhMNX1ozbsrmdtz58suYD1NGjR8kayNsG531wwKuLgw1NHOjq1gjwMjjw+Oabb/R+hu5Bva4aNQ6QG1t2bhPnnmP61JcbA9o0e2QxHs/GkMHSGsK1rTt27BB5Ebzv8HGKf7cxY8aotz15Pq5p4BrhjRs30ltvvSXyjzinjfM5TM2YvfjAOthF8FEf7o7JOyFfxckRNONByLg6XHfALK6y5L9XXnlFXLVz9S1fVXCyVV0BjHz1wFfO3ARkzLIzTjLVvELhKmlOOGzqZ2kuVxcPwsPJhsbo3slJk3wVzoP38IGuLnzw44MTD0xmSNLpddddJ5LJdu7cSf3792/y9+dtQrNGgpvbeL1yoGQschIiBw9N/b14GdwUwcnAxgoweZkNBXA8DycKjxgxoknV7HUtk/fH48eP1xloav4+mtt9ZmZmrS7aPK9cu6RJ37ZtbI1ZJ9z8pklO/K2Pvu/Fib5ygjKvi7///lvUfHAzVn3vkwMgTnTnP64l4URTPs7VFXzwb8UXYVy7Yozaj4bwb8nbOddqatZ+8DFJft0UjLVtgza7yPmoD1cvsnfffVdrunw1N378ePWOrHtlKB8cy8rKxKPc80B3aGU+sXAzzSeffCLyKXTxQbMp+GTFzSjvv/++Vtm4pwZXqcplbyw+CPF34x5Amt+FT0Z8VSSvM2Ps1Fx2bjOuL6jgK2jOY+HP5h43uvhkxTVQ3AuDPf744yI44oBQ30i2XF3MPQPqwrVkXFvAtTKavSc4q9/Yw2ZzTRjX1nAvCn25RIZsG3zVyr0FFi5cWOs17n3QlDJzkwsHFnJvL03ytsafy7lKPBKkLu4dpdlLi3MN+CTBPXnqw02ZXLPDvRA0r8w1P5e3ew7keVvQ3O5192HG2yrXTO7Zs0drndZVS2RMvA0auu75O2n+6daE1JUfo5krxt+R8zbkYEGuTdA9bumuJ952NJtg5GMW1yDKx7a6thFeNgc3za1FMwT/llxW7jmliZvB+Vhiqhqauo7r0Dx2X/PBV7Hc1slXyrxxcW4H78R84uUDIXfTY/z/jz76SCTgccTP0TcfdPnEIZ+M+aqTmw5++OEH0fbKVwM8TDj/cd4JX+lzNTWfSPmKjU+MfHXOJ00+0DcWnyCfeuopsfNzNSp3Q+MrOi4ndwXk7rRNxdWuvDNzzQF3a5O72nJbNDeTGAt3ueW/hnBwwUEDJ5txVTDXbnCb8Pnz50XiI5/YOJmY8e/DtVJTp04VtVncds2/AQcSXAXN89dXpc0nNu7Gyd0rueaDl8M1HlwVbuycD95+uFsfB198pcnfgX9X/l6cMMy1GboHW128zXJZuXs4JwpzMi9/B77C5e/KgRZ3rW0MrqbnZkbOy+HmRc4D4Ctc7mrLQRnvN1xmzi/gREFOLuWy8smBfwuezlfFcnOnoV1tufmRu0JyIMXNANy1kZuauOs4nwz5O/L64SZMfs7bAe9/3KWYEw41B62TA1Fu0uL9Y/78+equtnyVrNul29h4nfGVOl/IcNk5j0tf/klT8briYwonIHOQwEEFDx/A31netoYMGSLyOTiw5ZwvDuB5W9bExzLOI+JthH9XbprhcvM65/2uLnxs5G2ALyB4W5ObcrirLb/G3XaN6frrrxfL5e2DhzngsvL34cEiOdnY2F2ZZfUd16EZJBtWV3cqXTzI2IsvvigG4uGhtyMjI2sNMnbgwAHp1ltvFQP6yIM4XXfddWLQLk07duwQ3cO4K6Bu96zk5GTRTbVFixbic3jAI17GypUrDSpzXYOMcdfa9u3bi2XyIEI8OFldg4w1BndRHjhwoOjS5+vrK11//fVag4w1p6ttfeoqK3c9/vzzz8WAU9w1kL8vd/3jQZP0dcM9deqU6KrKg8Hx7+Hj4yO+z5IlSwwaJpoHQZMHZ+rVq5dBg4w1thuf5rrhwZn4e/GAUTxo18yZM7W2L3kwprpwV2De9vj34u/KXbsff/xx6dKlS7UGGdOl+70Yd7N+4IEHxHYqDyDGZdAcoIq7Jb/xxhvi9+L1xAPocRl4f+Juz43taivjbuLc1VJeJpdNcyhs7obLn8GD6jU0yNjhw4fF+5szyFhTutomJSWJgb+4fKYaZGzRokXieMXrifcL7g6vibvicrdl7s7N29aUKVPE9qB5bOLu5NzVPSEhQWw3vI3xc91BAPUNMsb7JJeDjz/y4GQ8AJ7m4GQNqWu7lrcZTdw1nbsSh4eHi/2fhylo7CBjTdlH6zuuQ9M48D/NCV4AAAAAGsPucz4AAADAvOw+5wMA7BsPJc9/9eE8E3QHbRxOYuVcsfro615ubZ8JTYNmFwCwa/L9kerDSZqa99iBhnFSNyfq18fYpx8lPhOaBsEHANg1HsJbcxhvfbhXCY/CCYbjsVoaGq7dmOMeKfWZ0DQIPgAAAMCskHAKAAAAZoXgAwAAAMwKwQcAAACYFYIPAAAAMCsEHwAAAGBWCD4AAADArBB8AAAAgFkh+AAAAACzQvABAAAAZoXgAwAAAMwKwQcAAACYFYIPAAAAMCsEHwAAAGBWCD4AAADArBB8AAAAgFkh+AAAAACzQvABAAAAZuVs3o8DsG/5pRX0074LtOVUJl3ILSZPVydq38KXBrcJpjGdW5Cbs5PSRQQAMDkHSZIk038MAKw9cpme/fUoZReV63090MuV7hvammYOjCEXJ1RKAoDtQvABYAYf/nOG3tpwUjyPC/Gi2/tGU4eWvqIm5FDaVfrlwEVKzy8Vr7cO8aIXJ3SmQW2CFS41AIBpIPgAMGPgcc+QOHp0dLtaNRuVVdUiAHljfZK6ZuSmHhH03HUdyc/TRZFyAwCYCoIPABP6/dAlevC7RPH8sdHtaO7w+HrnzyupoHf+PEVf7Uwl3jNDfdzoxQmdRD6Ig4ODmUoNAGBaCD4ATCQ5s5AmLNlGReVVdO/Q1vTk2PYGv3f/uRx6bOVhOptZJP7fI8qfHrgmnoa2DSUnRwQhAGDdEHwAmEBFVTVN/GA7Hb+cT31jA+mbO/uScyOTSEsrquijf87QZ1tTqKSiSkwL93OnER3CaGB8MPWKCaBgbzcTfQMAANNB8AFgAp/+m0yvrk2iAE8X2rBgCIX6ujd5WRn5pfTpv2dp5YELdLW4Quu1mCBP6hkdSNd1bUlD24aQI2pFAMAKIPgAMLKLV0vo2kVbRG3FW5O70pRekUZZLteEbDudRf+ezqRdZ7Pp1JVCrde5l8xjo9vT6E5hyA8BAIuG4APAyO7/Zj+tPZJOfWIC6Yd7+pksEMgrrqADabn076lMWrnvAhWUVYrpt/aJooUTOzW6mQcAwFwQfAAY0dGLeXTdkm3E8ca6+YPF6KXmUFBaQUs3J9PSLcmil8yk7q1o0ZQENMMAgEXCpRGAES3aqBrPY2JCuNkCD+bj7kKPj2lPS6f1JGdHB1qVeFEEIgAAlgjBB4CR7D+XS/+czBRdYRdc21aRMvB4IAtv6KwOhBLP5ypSDgCA+iD4ADASbvZgk3tEUEywl2Ll4JyPid3CqVoienzlYdHtFwDAkiD4ADCC1Kwi+jvpinh+99A4pYsjRkXlG9Wdziikb3efV7o4AABaEHwAGMGy7Ski0XN4uxBqHeKtdHHI39OVHhqpavr54J8zVFapGqQMAMASIPgAMEKX1x/3XRDP7xysfK2HbGqvSGrh606ZBWX0a+JFpYsDAKCG4AOgmb7be14MKNa+hQ8NaB1ElsLV2ZHmDIoVzz/59yxVcxIIAIAFQPAB0AyczPnVjlTxfPagWIsbWfSWPpHk4+4sblD3d1KG0sUBABAQfAA0wz9JGXQ5r5SCvFxpQkI4WRoe/+O2vlHi+be7zyldHAAAAcEHQDP8uC9NPE7uGUHuLk5kiTj3g205lSluUgcAoDQEHwBNxCdyHlSMGevmcaYQF+JNPaL8xbgfvx5E4ikAKA/BB0AT8S3uq6ol6hUdQPGhynevrc/knqrgaOX+C4TbOQGA0hB8ADQBn8B/qulee7MF13rIxndtKXq/nLpSSEcv5itdHACwcwg+AJpgT0oOpWQVkZerkzixWzo/Dxca2SFMPP/jyGWliwMAdg7BB0ATyIOKXdc1nLzcnMka8E3n2Pqjl9H0AgCKQvAB0Egl5VXiBM5u7h1B1mJ4+1DR9JKaXUwnrxQoXRwAsGMIPgAaaVNSBhWVV1FEgAf1iAoga+Ht5kxD2gSL5+uOpCtdHACwYwg+ABrpt5ruqjyomKWNaNqQMZ1V+SnrjyL4AADlIPgAaIS8kgraXDO2x8RurcjacNKps6ODaHY5n12sdHEAwE4h+ABohA1H06m8qprahflQuxY+ZG38PF2oR7SqqWjzKdzrBQCUgeADoBF+P3RJPE7oZnn3cTHUsHYh4lGuwQEAMDcEHwAGyigopR3JWeK5Jd5EzlDD2oaKR/4upRVVShcHAOwQgg8AA609fFncH6V7lD9FBnqSterQ0ofCfN2otKJaDJYGAGBuCD4ADLSupofI+C6WP6JpfbiHztC2aHoBAOUg+AAwQFZhGe1NzdEaKdSaDWunanpB0ikAKAHBB4ABNh67Ippcukb4UUSA9Ta5yAbGB5OTowOdzSyitBx0uQUA80LwAWCAdTXDqdtCrYd8o7mecpfbk6j9AADzQvAB0IC84gramZwtno+tGSHUFqDLLQAoBcEHQAP+PHGFKqslat/Ch2KDvchWDGmjCj52nc2miqpqpYsDAHYEwQdAA+Q72I7uZHlNLseOHaMpU6ZQSEgIeXh4UJcuXejdd9+l6uqGg4mOLX0pyMtV3CQv8fxVMW379u00btw4CgwMJG9vb+rTpw/93//9X73LWbduHY0cOZL8/f3J09NTlOHNN9+kioqKOt9z4cIFuvfeeykqKorc3NwoPDycZs6cSSkpKU1YCwBgbRwkSZKULgSApSoqq6TuC/+k8spqWr9gMLVv4UuWYufOnTRixAgqKSkRQUJMTAz9+++/lJ6eLgKSH374ocEb3837LpFWH7pE866Jp5iCozR16lQRuAwZMoSCg4Pp77//pqtXr9IjjzxCb7/9dq33v/HGG/Tkk0+So6Mj9e3bl0JDQ0W5MjIy6Nprr6W1a9eSi4uL1nuOHj1Kw4cPp6ysLFHmHj16UHJyMh06dIh8fX3Fd0hISDD6+gIAC8LBBwDot+7IJSn6iTXSkDc3SdXV1ZKlKC8vl2JjY/nCQVq8eLF6ekFBgdS/f38xfdmyZQ0u54c958X3G/vGWsnX11e87+eff1a/np6eLsXHx4vp//zzj9Z79+zZIzk4OEguLi7S+vXr1dPz8vKk4cOHi/e8+uqrWu/hddilSxfx2uzZs6WKigr1a++//76Y3rFjR6mysrLJ6wYALB+CD4B6PPzDQXFyfmn1McmS/PDDD+JEnZCQUOu1/fv3i9c6d+7c4HIu5haL7xcwbKZ4z8SJE2vN88svv4jXrrvuOq3pc+bMEdPvuuuuWu85efKkCEyCg4O1AomtW7eK9wQGBopASdeAAQPE67/++muDZQcA64WcD7AIqampoolg2LBhVFRURA8//DBFRkaKPAaull+9erV63p9++klU8Xt5eVFYWBg9+OCDoulBV3FxMb322mvUvXt3kb/Af/369aOvvvpKbxm2bt1KDzzwAHXt2pUCAgLEZ3849zrK3byc+rZyrzX/5s2bRZk5VyEnJ4fuu+8+atmypchh6Ny5M3355ZdkKn/88Yd4nDx5cq3XeH3FxcWJ5g1er/UJ9/eg1iFeVJy8r87ljR8/ntzd3emvv/6i0tJS9fT9+/eLR/7NdLVt21bkcXDTCueR6L6nZ8+e4vfQxc0x7Lfffqu33ABg3RB8gEUpLy8XeQzffPONCBT4j3MBJk2aJE5+77zzDt12223k4+NDo0ePpqqqKlqyZAndeeedWsvhnIP+/fvT008/LXIghg4dKvIYkpKSRLAwb968Wp/92GOP0RdffCGCDi5DzwFDqLK0mPJ3r6SH7phAhYWFesvMORH8Wb///jsNHjyYBg4cKD5nzpw59Pnnn5tkPfE6kQMNfeTphw8fbnBZg9uEUHlGSp3Lc3V1FcEUBx6nTp1ST+cgkXGgpk9QUJBWWZv6HgCwQUpXvQCwlJQUUd3Of9dcc41UWFiofo1zF3g65x4EBARIe/fuVb928eJFKTQ0VLyenJysnj5u3Dgxbf78+VJpaalWDkOvXr3Ea+vWrdMqw9q1a6WrV6+q///qH8elqEdWSZ1H3Cjmf/HFF7Xm5xwIucy33HKL1uesWrVKTI+Kiqr1XYcOHap+n6F/uvkbvB54+qFDh/SuzwULFojXOY+iIb/uPq3+HM7X0OeGG24Qr//++++1mkiWLl1aa37O7fDz8xOvP/LII+rpn376qZjWt29fvZ8zb9488XpQUFCD5QYA6+WsdPADoIl7TSxdulQ0qcimT58uaiXOnDlD//vf/6hXr17q17hqf9q0aaJGhHtJcHPDwYMHRS+L3r170+LFi8UyZdxM8+mnn4orfP6cMWPGqF8bO3ZsrfE9HJxd6IVX36RbBv4umgKee+65WmXmHhoffPCBaG6R3XDDDaK2QG764F4dMv5Mzf8bIj4+Xuv/ci0Md23VR15/BQUFDS67Y4ir+nlWCX8fw5bHNUk7duwQzVjcbVbTzz//THl5eXrfw/bu3UvHjx+njh07ajWT/fjjjwaXGwCsF4IPsCh8UuZ8AU0cPERHR4v8gVGjRtV6Dwcc7PJl1XgcGzduVAcAmoGHTM4B2bNnT63XLl68KPJLdiceoT07T5EjSfTLxXDR9HD69Gm9Zeb8Bbm5QBN/Dw4+uFyawQZ3TbUkXm7/HQa2ncmkuDDDuhPff//99OGHH9KuXbtEgPjss8+K8UZ4/XP+i7OzM1VWVmr9Bu3atRNNaKtWraIJEyaIQJCDxLNnz9KCBQsoO1s1kqy+3w0AbAeCD7AorVq10jtdTk7U97r8WllZmXiUkyyfeeYZ8VcXzeRJxrUkHBjoDo717ZH6yxwREaF3OuelaJbLmPg75+bmitoCfeTcCrkMDS1LtuX4BZo+sLVBy+OE4F9++UWMKbJixQrxJ+OAgoM8Di508zs4r4aDDK6p4twaGS+bByfjZOO6ckIAwDYg+ACL0tAVryFXxPLonoMGDaLWrWufSPXhq3ceSMvPz4/ee+89+umiDx3JcaQXJyXQzIGxonlHrllpSpk0vf766yIhtTE4oZa/j4xHBuXgg0cK5d45ung64xqjhnCzkbevLxXm59O2g6eoqnqIuOOtIcvjgcS41uL7778XtTxOTk40YMAAuummm2jWrFlink6dOmm9hwML7inEI6PyIzfP8O/EzWcnTpzQ+x4AsC0IPsDmyDUR3OzCAYUhuBmAvfLKK3Tj1Nto4Ut/koOzRNe0DxPdeLnHjLGsX7+etmzZ0qj3cHdWzeCDRwDlHiEHDhwQw6Hr4ulMX2CiT/eEbrR167+Uff4kHbmYR90i/dWvcU0QBxbc3Va3SUwOJriZRRePdMqBmZznoYm7KHO5dcvOtSLy9wUA24WGVbA5fJ8RzYDCEFyLIAcuO5KzxY3k+CZyUUGeYlwRY96FgK/2awb4M/iPuwfrjr3BVq5cWWv5iYmJojaCE14NTWy97jrV8opPbqdtp7XvcrtmzRrRRMW1HByAGDoOCZeBk2u5ecYQ3ITEwQfn18yYMcOg9wCAdULwATaHByDjAIQHt5o7dy7l5+fXmodrDbgGQiZf0fPJb9PxS+L50LYhokfGE088QZaGkzZjY2PF9+CePpq5Gfydmb5aH86xaN++fa1kW27W8fDyoZLTu+i7H1dqjZfy+OOP17k8HjRMNzDjHjDc5MKBCufR6OKxQnR/Ex6kje8rc/78eTE2S115NABgG9DsAjbp66+/FlfdH330EX377bfUrVs3kbfB+QU88FZaWhrNnz9f3dWWT5aLFi0SPV02bttHTmHxtOZfB3p1307RfMMn63PnzpGl4Ju18Xfk2ghO0OSbyHE+Bo/SyrkpPFKpvtoDvoEbfw/dRFW+i+2iJUvp/jnTaccnz9DQpI0UGhIsBnbjQdT4M/Q1hXBuBw/0xrUs3PzCPYI4IOHAg2tluHeLLv49+IZ0nJTKCcT8m3C5ufsw1/BwrxkAsG2o+QCbxHdX5Svw999/X4wlwU0RfDLkwIO75r711lv06KOPqufnrrI89sT1N95MlRXlVHJmNxXlZtLChQvpu+++I0vEiZ1cZg4AeAwUHmGVgwiubTDkjra67p15G3W55x1yj+1BBxITxVgpPL7I8uXLRWCm9z333isCCA7OeGyPzMxMuuuuu0SOiNw0pOuaa64Ro9OmpKSI3jK7d+8W34WfL1u2DN1sAeyAA480pnQhACzFF9tSaOGa4zS4TTCtmNOX7M1Tvxyh7/acp1kDY+j569HjBABMA5cYABq2nMpU53vYIw662LbTWUoXBQBsGIIPgBol5VW062y2XQcfA1oHEbfWnM4opPQ87UHYAACMBcEHQI1dKdlUXllN4X7uFB9a+3bv9sDf05W6tvITz7edQe0HAJgGgg+AGjtqTrZ8i/nGJmvakkE1TS9bdcb7AAAwFgQfADV21jS5DIivfZM4e8LBF9t+Jouqq5GPDgDGh+ADgIiuFpfTsUuqga/6x9l38NEjKoA8XZ0oq7CcktJxa3sAMD4EHwCc73E2h7jTOed6hPoaNoS4rXJ1dqS+sYHi+bYzaHoBAOND8AHATS7JWereHsB5H6qml63ocgsAJoDgA0Aj38Pem1x0x/vYk5JDpRVVShcHAGwMgg+we5kFZXTqSqF43g/Bh9Am1JvCfN2orLKa9qWq7vgLAGAsCD7A7sm1Hh1b+lKAl6vSxbEI3NV4UHxN0wvyPgDAyBB8gN1Dvod+GGodAEwFwQfYvZ3JNfkeCD60DIxXBR/cBTm7sEzp4gCADUHwAXaN71+Sml1Mjg5EvWu6l4JKiI8btW/hI55vrwnQAACMAcEH2LW9qTnisUNLX/J1d1G6OBbc9IK8DwAwHgQfYNf21QQfvWNQ69HQeB8Sj8IGAGAECD7Aru2t6UaK4EO/PjGBYsTTy3mllJxZpHRxAMBGIPgAu5VfWkEn0lX3c+kVE6B0cSySh6sT9a5ZN2h6AQBjQfABduvAuVxxP5eoQE8Ks/P7udRHHu9j2xl0uQUA40DwAXZLHrkTTS6GJZ3yzfcqqqqVLg4A2AAEH0D23tNFblYA/Xjk10AvVyosq6SDaVeVLg4A2AAEH2CXyiqr1CfSXqj5qJejo4N69Netp5D3AQDNh+AD7NLRi/nipml8Rd86xEvp4li8IXKXW+R9AIARIPgAux7fo1d0gLiJGtRvUE3ex6G0q5RXUqF0cQDAyiH4ALu075wq2RRdbA0T7u9BcSFeVC39dy8cAICmQvABdodH6kw8r8r36BmN4MNQg2tuNLftDPI+AKB5EHyA3bmQW0JZhWXk4uRAncL9lC6O1Q21vu008j4AoHkQfIDdOXBe1eTSMdyP3F2clC6O1egXF0hOjg7iLsBpOcVKFwcArBiCD7A7cpNL90h/pYtiVXzcXahHlL/6RnMAAE2F4APsTmLN+B7da06k0JSh1pH3AQBNh+AD7EppRRUdv5QnnveIQrJpU7vcbj+TTVXc9QUAoAkQfIBdOXYpnyqqJAr2dqWIAA+li2N1EiL8yMfdWYz1cfSiKogDAGgsBB9gVxJrkk27RWJwsaZwdnKk/nGqodZxl1sAaCoEH2CfyabI92j2XW63nkbeBwA0DYIPsMuaD+R7NN3gmvE+9p/LpaKySqWLAwBWCMEH2I30vFK6lFdKjg5EXSMwuFhTRQd5inwZzp3Zk6K6Rw4AQGMg+AC7cTBNVevRroUvebk5K10cq8W5Mv81vSDvAwAaD8EH2A3kexjPwJr7vOw6i5vMAUDjIfgAu4GRTY2nd0ygeExKz6eC0gqliwMAVgbBB9iFiqpqOnxRrvlAsmlzhfm6U2SgB/E4Y3JQBwBgKAQfYBeSLhdQaUU1+Xm4UFywl9LFsQm9o1W1H/tSkXQKAI2D4APsQmJNsmm3SH9y5O4u0Gy9appe9p1TrVsAAEMh+AC7gGRT4+sVE6Bet9ysBQBgKAQfYFeDiyHfw3jiQ7xFM1aJuFlfvtLFAQArguADbF5OUTmlZheL590iUPNhLNx81TNaFcyh6QUAGgPBB9jN4GKtQ7zIz9NF6eLYZNMLkk4BoDEQfIAd5XugycXYesk9Xs7lkiRJShcHAKwEgg+weUg2NR2+R46rkyNlFpTR+RxV0xYAQEMQfIBNq6qW6GCaPLIpaj6Mzd3FiTq38hXP96Ui7wMADIPgA2xacmYhFZZVkqerE7UN81a6ODapW01Qd/gCRjoFAMMg+AC76GLLzQPOTtjcTSEh0k88HrqQp3RRAMBK4GgMNg3JpqbXtab78vHL+VReicHGAKBhCD7ApuFOtqYXE+RJvu7OIvA4daVA6eIAgBVA8AE2i3M9TmWoTobd0NPFZBwcHCihJrg7hLwPADAAgg+wWYfTrhIPPRER4EGhPu5KF8emcU4NO5yGvA8AaBiCD7BZiTVdbPlOtmCevA/UfACAIRB8gM3CzeTMJ6Em+DidUUgl5VVKFwcALByCD7BJPNQ3RjY1nxZ+7hTq4yYGdTt2CU0vAFA/BB9gk9JySii7qFwM/d0pXDUCJ5ir6QXBBwDUD8EH2KTEmjvZdgz3JTdnJ6WLYxcS5KRT5H0AQAMQfIBNQpOL+XWtSew9jJoPAGgAgg+wSUg2Nb+urVQ1HylZRZRfWqF0cQDAgiH4AJtTWlFFxy7li+cY2dR8ArxcqZW/h3h+omb9AwDog+ADbA73tqislijY200MMAbm06Glr/o+LwAAdUHwATad78FDf4P5cIIvO46aDwCoB4IPsDn7UlX5Hj2Q72F2crdmudkLAEAfBB9gc4OL7TunCj56xSD4MLeONc0upzMKxF1uAQD0QfABNuV8TjFlFZaJwcW61PS+APPhHBsfd2eqqJLoTEah0sUBAAuF4ANsssmlcytfcnfB4GLmxjk2cu0Hkk4BoC4IPsCm/NfkEqh0UexWp3BVjRPu8QIAdUHwATZl/7kc8dgrGvkeSkGPFwBoCIIPsBlXi8vp1BVVnkFPBB+K0Wx24QRgAABdCD7AZhyoGVI9LtiLgrzdlC6O3YoP9SYXJwcqKK2kC7klShcHACwQgg+wuWRT1Hooy9XZkdqG+YjnGO8DAPRB8AE2A+N7WA70eAGA+iD4AJvAA1odSlMNq94zGj1dlIakUwCoD4IPsAncrbOsspoCPF2odYiX0sWxe+1aqJpdTl5B8AEAtSH4AJuw/9x/+R64mZzy2rdQ1Xyk5ZRQYVml0sUBAAuD4ANswt7UmvE9MLiYRQj0cqVQH1WPo5PpBUoXBwAsDIIPsHo8loRc84HBxSyw6QXBBwDoQPABVu9cNt9MrlzcTK4zbiZnMdqrgw/kfQCANgQfYDNdbLtE+OFmchakXU3eRxJqPgBAB4IPsHq4n4uF13xcKcAw6wCgBcEHWD2MbGq5w6w7OvA9dyooo6BM6eIAgAVB8AFWLaeonE5nqG4mh54uloWbwGKDVWOuoOkFADQh+ACb6GLbNsxbdO8EyxzvA0mnAKAJwQdYtT0pquCjN2o9LLq7bdJl1HwAwH8QfIBNBB99YhF8WHTwgWYXANCA4AOsVkFphbinC0PwYdk9Xs5kFlJlVbXSxQEAC4HgA6wWj2paLRFFBXpSSz8PpYsDekQGeJKnq5O463BqdpHSxQEAC4HgA6wWmlwsn6OjA7UNQ9MLAGhD8AFWC8GHtQ2zjuADAFQQfIBVKq2oosMXVPkefRF8WDQknQKALgQfYJUOpl2l8qpqCvN1EzkfYA3BB8b6AAAVBB9g9eN7ODg4KF0cMGCgsbScEiosq1S6OABgARB8gFUHH2hysXw88myoj5t4fuoKml4AAMEHWKGKqmrRzZb1iQ1SujjQiKYXJJ0CAEPwAVbn6MU8KqmoIn9PF2oT6q10ccAA6PECAJoQfIBV53vwOBJg+drV5H0g6RQAGIIPsDrI97Dumg9JkpQuDgAoDMEHWJXqaon2pmJwMWsTH+pNXEmVW1xBGQVlShcHABSG4AOsyskrBZRfWklerk7UsaWqKh8sn7uLE8WFqPJzMNgYACD4AKtscukZE0jOTth8rbHp5cRl5H0A2DscvcE67+cSE6B0UaCROtTUVCUh+ACwewg+wGpwouJu9c3kML6HtdZ8oNkFABB8gNVIySqirMIycnV2pK4RfkoXBxqpfU3Nx5mMQiqvrFa6OACgIAQfYHVNLt0i/UUCI1iXcD938nF3pspqiZIzC5UuDgAoCMEHWA2M72Hd+AaAHTDYGAAg+ABr8l++B4IPa9W+ZU3ex2XkfQDYMwQfYBUuXi0Rf06ODtQjCj1drFX7mpqPE0g6BbBrCD7AKuxKzhaPXVr5kZebs9LFgWbXfKDZBcCeIfgAq7DzrCr46N8aXWytWbswVfDBQ6xnF2KYdQB7heADrMKumuCjXxyCD2vGtVbRQZ7qm8wBgH1C8AEWLy2nmC7kqvI9ekUj38NmhllH8AFgtxB8gNX0cuGBxZDvYTtJp8j7ALBfCD7A4u2sSTZFk4tt6CAnnaLmA8BuIfgAq8n36I/gw6ZqPk5dKaDKKgyzDmCPEHyAxed78Pgezo4O1BP5HjYhKtCTPFycqKyymlKzi5UuDgAoAMEHWEUXW+R72A5HRwdqp77DLfI+AOwRgg+waOhia9t5HyeQdApglxB8gMWSJIl2n1X1dMHgYralQ0tV3sfxSwg+AOwRgg+wWDy2B/I9bFOncD/xeBTBB4BdQvABFt/FNiHSnzxdke9ha80ujg5EmQVllJFfqnRxAMDMEHyAFeR7BCpdFDAyDiZbh3iL50cu5ildHAAwMwQfYLH5Hkg2tW18h2J29CKaXgDsDYIPsEhpOSV0Ka+UXJyQ72GrOsnBxyXUfADYGwQfYJF2ns0SjwkRyPewVZ3DVT1ejqHZBcDuIPgAi7QD93OxeR1rgg+u4couLFO6OABgRgg+wCLzPbafUdV8DGoTrHRxwER83F0oNthLPD+GLrcAdgXBB1ick1cKKKuwnNxdHKl7lL/SxQET6oy8DwC7hOADLM72M6omlz6xQeTm7KR0ccAMeR9HkfcBYFcQfIDFUTe5xCPfw25qPtDdFsCuIPgAi1JRVU27a8b3GNAa+R62rlNNzcf5nGLKK65QujgAYCYIPsCiHEq7SkXlVRTg6UIda24+BrbL39OVIgI8xHPkfQDYDwQfYFG21TS5DIgPJke++QfYvG6RqqTixPO5ShcFAMwEwQdYlB01yaaD4tHkYi96RKlGsE08f1XpogCAmSD4AItRVFZJB2qufgci38NuyN2pE9OuijFeAMD2IfgAi7EnJYcqqyWKDPSgqCBPpYsDZhzp1NXJkXKKyulcdrHSxQEAM0DwARbYxRa1HvaEx3Lp3EqVXCzXfAGAbUPwAZaXbIomF7vTHXkfAHYFwQdYhIyCUkpKLxDPB7TG4GL2mnSKmg8A+4DgAyzCv6dUtR5dI/woyNtN6eKAQkmnHIAWl1cqXRwAMDEEH2ARtpzKFI9D24YoXRRQQLi/B7XwdaeqaokOX8BgYwC2DsEHKI5POFtPI/iwd+out8j7ALB5CD5AcYcuXKWrxRXk6+6sHu0S7Df4QN4HgO1D8AGK23JSVesxuE0IOTthk7RX/410movBxgBsHI70oDjkewDrEuFHbs6OlFVYTsmZRUoXBwBMCMEHKIpHteRmFzYEwQfZ+2BjctPLrrOqe/wAgG1C8AGK4kRTrmFv38KHWvi5K10cUFi/ONUYL7tTcpQuCgCYEIIPUBSaXEBT39ia4ONsNvI+AGwYgg9QTHW1RP/KwUc7BB+g6vHCN5nLKCijVNxkDsBmIfgAxRy/nC+SCz1dnahXdKDSxQEL4O7ipO5uzbUfAGCbEHyAYjafzFDfSM7VGZsiqPSLUwWiOxF8ANgsHPFBMX8nqYKPYWhyAQ0D4lV3Nd5+Jks0zQGA7UHwAYrIKiyjg2mqLrYjOoQqXRywsMHGuCmOm+ROpOcrXRwAMAEEH6CITUkZoottp3BfaunnoXRxwIJwE1z/mi63W0+r7nYMALYFwQco4u8TV8TjiA5hShcFLNDgNqqmF/mGgwBgWxB8gNmVVlSpr2ivRZML6DG4ZtyXvSm5VFxeqXRxAMDIEHyA2fHolcXlVRTq40adw/2ULg5YoLhgL2rl70HlVdW0+yxGOwWwNQg+QMEml1BydHQwyWecOHGCpk2bRi1btiQ3NzeKiYmhBx54gLKy6s8hWL16NQ0dOpR8fX3F37Bhw+iPP/7QO291dTU999xzFB4eTh4eHmLew4cP6523srKSunTpQgMGDGjSyJ0ODg7irz7Lly8X88ycOVPvdM0/Ly8vUW4u8xNPPEHHjh1r9HJNiT9vSNtgrS7ZAGA7EHyAWfGJ9+8TqpPJiPamyffYtGkT9erVi7799lvy9/en6667TgQgH374IXXv3p0uXLig933vvvsuTZgwgXbs2EEDBw6ka665hvbs2SPe/8EHH9Sa/4033qCFCxeSn58fjRw5knbu3EnXXnstFRQU1Jp3yZIldPz4cVGGhoIIU2ndujXNmDFD/E2cOJE6d+4sgo4333xTPL/99tspP99yepdcU7N9/HWCk5PR5RbApkgAZnT8Up4U/cQaqe0za6XiskqjL7+oqEgKCwvjM5X03HPPqadXV1dLjz76qJg+atSoWu9LSkqSnJycJDc3N2nHjh3q6SdPnpSCgoIkZ2dn6fTp0+rp5eXlkr+/v5SQkCCVlpaKaV9//bVY/ltvvaW17PT0dMnX11e67777mvy9eLkN7a7Lli0T88yYMcOg6fJ6Wb16tRQTEyPmGTp0qPhuhr7flHj74O2Et5cTl/PM+tkAYFqo+QBFmlwGxQeTh6uT0Zf/yy+/0JUrV6hdu3b0/PPPq6dzbcOrr74qml82btxIhw4d0nrfe++9R1VVVXTvvfdS//791dPbtm1LzzzzjGg24XlkqampdPXqVbrllltErQq79dZbyd3dnQ4ePKi17Mcff5xcXFzo5ZdfJkvD64Vrdnbv3i2aYbZs2UJLly4lS8Dbx8CaAcfk2jIAsA0IPsCsuArdlF1s9+/fLx6HDBlCjo7amzcHANycwn777Tet1+S8jsmTJ9dapjyN80Fkubm54jEgIEA9jT+Pm2Dk1xg34axYsYJee+01Cgy03PvXhIaG0ksvvSSev//++2Qp5AHo5KAVAGwDgg8wm8yCMjp0wbSjmhYVFdUKCjQFBakGr9Ks+eAajPPnz4vnnBOiKzIykoKDg+ncuXPqnIioqCjxeOrUKfV8HHRkZmaqX+OEVE5y7dmzJ82ZM4cs3c033ywCqOTk5DrzYsxNzgtKTLsqRsUFANuA4APM5q8TV8SopgkRfhTm626SzwgJUY0PwYGCPikpKbVelwMPDli4F4g+ERERWu9r0aIF9ejRg5YtW0bbtm0TgcfDDz8sAo7x48eLeT7++GPRBMNJprq1MJbIx8eH4uLixHNOjrUELfzcqUsrP7HdbDyG2g8AW2H5R0SwGX8eV508RnVqYbLP4OYWuRlFt1vtxYsX6c8//xTPNXukFBYWikdPT886lysHJZrvW7RokahpGTx4sGhS4S6p48aNEzkU2dnZ9Oyzz9Ls2bOpT58+6veUlpaKAKWpdLvMav7NmjWLmotreJhm05HSxndtKR7/OHJJ6aIAgJE4G2tBAPUpLKukbWdUwcDIjqYbUn3UqFGiRuLAgQM0duxYUevQsWNHOnLkCN1zzz0icZQZoyaCx8jgz+GcDm666du3L91xxx3itaeeekp0D3399dfF///++2968MEHRY0CjwnC83ECKyeoNgZ3k63LmTNnaPv27c36TnKXVqW6A+szvktLen1dEu1MzhZNL8HeqgRfALBeCD7ALP49lUnlldUUE+RJbUK9TfY5fNLkHi/c9LFv3z4REMjCwsLohRdeoP/9739aOSHe3qryFBcXN5hLwk0Tmjp16qQOMGT8uV988YVI3OSaBK5xuf7668VYGj///LMIQLgcXJuyePHiRn0/rl2p77XmBh9ybZElJcdGBnqKprpDF/Jo/dF0ur1ftNJFAoBmQvABZrHxWLq6ycXUV9XR0dEi12LVqlWit0lJSYkIEnjEUw5MGP9fJieIclMDBxn68j7kBExedkM1B3PnzqWuXbuKbruMa1+4ueXHH38UXX1vvPFGUUvB07n7bX3NPebEybRnz54Vz7m2yJJw0wsHH38cvozgA8AGIOcDTK6iqpo2Jam62I4yYZOLJmdnZ5oyZQq98847IvFz3rx54mqegxG5yUTGo6DKAUhiYmKtZaWlpYkaAQ48eMj1+nz55Ze0d+9eMSKqk5NqHJOkpCRRA8KBh4zzQMrLy0UQYik4OOLgicc24TE/LMm4Lqq8j90p2ZRRUKp0cQCgmRB8gMntScmh/NJKCvJype5R+rvAmkN6ejqtXLlSdLfl2gdNcg8Vfl2XPI2bTurDeR+c68H5HPJ4IjKufdHXjGMpvWAyMjLEfWrY/PnzydJEBHhS9yh/qpaINhxV1aIBgPWyjCMf2EWTy7UdwsjJRDeS03T06FHRzKHbbML3M+HeKtxLhZM+NfEJl2squJZk165d6umnT5+mV155RdSkNHRS5lySsrIyca8UTdzEwz1q5IHNKioq6KeffhIjo/L9VpTENR1r164VuTGXL18W97O5++67yRJx4ilbfeiy0kUBgGZCzgeY/OT2Xxdb8zS5vP322yLfg3u98F1t+aqex+LgwIC7v+rrMcLDsb/11ltirA7uOss3inN1dRVDsXOtBSePxsfH1/mZPGgZBy782ZzYqolzQPimdVOnTqXRo0eLphZOOn3yySdrBUGmxOtAvjMtN/lwd2DurSMnmXKNDeehcKBliTjv45W1J2hPag6l5RSLRFQAsE6WeZQBm3HsUj5dyislT437dJjaDTfcIJpYOCDg3h/cs2XMmDG0YMECrVwPXQ899JAIMDgI2bp1q5jGd8fle7Pw2B314ZySDh06iBFNdfGAZBs2bKBHH32U1q9fL3JM+Lk8nLm58Mil/Mc46OFycGJpv379aPr06VpJuJaopZ+HuCfQ1tNZ9POBC7Tg2rZKFwkAmsiB7y7X1DcDNGTxxpP0/qYzNLZzC1p6e0+liwNW7tfEi7Tgh4MUFehJWx4bZlHjkQCA4ZDzASa1sabJxZQDi4H9GN2pBXm7OdP5nGLam2o5o7ACQOMg+ACTOZddREnpBSLJ9Jr2prmRHNgXD1cnGtdFNTz/yv1pShcHAJoIwQeYjJxo2jc2kPw9XZUuDtiIyT0jxePaI+lUUl6ldHEAoAkQfIDJm1zMNbAY2IfeMQEi54PvF7T+GLrdAlgjBB9gEtmFZbQvNUc8H2nCu9iC/eEk08k9I8Tz73aj6QXAGiH4AJP4OylDjEbZKdyXWvmbbywLsA8394oUuUQ85sfpKwVKFwcAGgnBB5jExmNykwtqPcD4Wvi504iaJOZvdp9XujgA0EgIPsDoOAlw25lMs45qCvZnWs3dbXnAMSSeAlgXBB9gdP+ezqTSimqKDPSg9i18lC4O2KjB8cFiGysoraQ1hy8pXRwAaAQEH2CyJpeRHVpgBEowGUdHB7q1T5R4jqYXAOuC4AOMqrKqmv5OMu+N5MB+TekZSS5ODnQw7SoduZCndHEAwEAIPsCouPfB1eIKCvB0oV7RAUoXB2xciI8bjevSUjxftiNF6eIAgIEQfIBRbTiarr6Xi7MTNi8wvVkDY8XjmkOXKbOgTOniAIABcHYAo6mulmj9MVXwMaYzutiCeXSL9KfuUf5UXlVN3+w+p3RxAMAACD7AaA5euEpX8svEXUcHxgcrXRyww9qPr3edp/LKaqWLAwANQPABRm9y4TvYujk7KV0csCNjO7egMF83yiosoz+OoNstgKVD8AFGIUkSrasJPtDkAubm4uRId9QMOrZse6rYHgHAciH4AKM4cbmAzucUk5uzIw1rF6J0ccAO8Zgfrs6OdPhCHh04n6t0cQCgHgg+wCjWH1Xd2nxo2xDydHVWujhgh4K83WhiQrh4/sU2dLsFsGQIPsAo0MsFLMGcwarE0/VH0yk1q0jp4gBAHRB8QLMlZxbSqSuF5OzoQCM6YFRTUE77Fr6i2a9aIvps61mliwMAdUDwAc3GV5lsQHww+Xm4KF0csHP3Dm0tHn/afwGDjgFYKAQf0GwbappcuLsjgNL6xgZSQqS/GO/j/3amKl0cANADwQc0y4XcYtG7gG9ey0OqAyiN76R839A48fz/dp6jorJKpYsEADoQfECzbDimuoNt75hACvZ2U7o4AMLIji0oNtiL8koq6Pu9aUoXBwB0IPiAZll9SDWaJJpcwJI4OTrQXYNVtR9fbD1LFVUYch3AkiD4gCZLyymmg2lXydGBaHxX1W3NASzFjT1aidq4S3ml6iAZACwDgg9ost9rDuj9WwdRqI+70sUB0OLu4kSzBsaI559sOYsh1wEsCIIPaDL5anJCzaiSAJbm9n7R5OXqRCevFNDmk5lKFwcAaiD4gCY5daWAktILyMXJgcZ0QpMLWCYed+a2vlHi+dItyUoXBwBqIPiAZtV68L1c/DwxsBhYrtmDYsXou3tSckSOEgAoD8EHNBq3ncv5HtejyQUsXEs/D5rQTbWdfvovaj8ALAGCD2g0vno8l11MHi5OGFgMrMLdQ+LUtwI4l40bzgEoDcEHNBrfM4ON7hRGnq7OShcHoFE3nPt8a4rSxQGwewg+oFFKyqto9UFVk8vNvSOVLg5Ao2s/ftqfRjlF5UoXB8CuIfiARll/7DIVlFVSZKAH9YsNUro4AAbrHxdEXVr5UWkFbjgHoDQEH9AoP+1TNblM7hFJjjy0KYCV4BvOybUffMM5rsUDAGUg+ACDJWcW0o7kbHEH25t6tlK6OACNxvcgigjwEM0uKw+oAmkAMD8EH2CwFTvPiccR7UMpIsBT6eIANJqzkyPdOShWPP9861mq4gxUADA7BB9gkKKySvq5ppfL9P6q+2UAWCNOlPb3dBHdxTccS1e6OAB2CcEHGGRV4kWRaBob7EWD4oOVLg5Ak3H38On9osXzT/7FDecAlIDgAxrEVdNfbk9R36gLiaZg7aYPiCE3Z0c6lHZVDLsOAOaF4AMatO7oZTqbWSRu0jUVY3uADQj2dqObekaI55/+e1bp4gDYHQQfUK/qaok+2HRGPJ89MJa83TCiKdiGuwbHiZ5bfydl0JmMAqWLA2BXEHxAvf46cYWS0gtE0DFzABJNwXZw/tLIDqp7E32xDUOuA5gTgg+oU2VVNb254aR4Pr1/NPl5uihdJACjuqtm0LGfD1ykrMIypYsDYDcQfECdftx3gc5kFFKApwvdM7S10sUBMLpe0QGUEOlP5ZXV9PUu1Tg2AGB6CD5Ar8KySlr85ynx/MERbUSyKYAtDrl+1+BY9SB6pRUYch3AHBB8gF5vbzgpqqGjgzxpWl/VmAgAtmhMpxbUyt+DsovKxXg2AGB6CD6glsTzufRVzV0/F07sTK7O2EzAtodcnzUwRj3kOvfwAgDTwlkFtHC185M/HyEe9PHG7q1oSNsQpYsEYHI8fo2PmzMlZxbRllOZShcHwOYh+AAtb64/SSevFFCQlyv977qOShcHwCx83F3olj6qAfQ+24pBxwBMDcEHqG0+maEeRv2tKV0p0MtV6SIBmM3MgbHk5OhAO5Kz6dilPKWLA2DTEHyAkJZTTA/9cFA8n9E/mq5prxp8CcBecNLp+C4txfPPt2LQMQBTQvABVFxeSXf93z7KLa6gLq386KlxHZQuEoAi7qzpdrv60CW6nFeidHEAbBaCDzvHtxN/7KfDYgj1YG9X+uSOnuTu4qR0sQAU0TXCn/rEBlJltURf7cCgYwCmguDDzn20OZn+OHKZXJwcaOntPSnc30PpIgEofsM59u3uc1RUVql0cQBsEoIPO7Yp6Qq9vVF175YXJ3Sm3jGBShcJQHEj2oeKm87ll1bSj/vSlC4OgE1C8GGn+J4t8787KMbzmNY3im7rG6V0kQAsgqOjA80epMr94N5ffINFADAuBB92KL+0gu5esY8Kyiqpd0wAPX99J6WLBGBRJveIEDdUTMspEc2SAGBcCD7sDA8dveD7g3Q2s4jC/dzpo2k9MXw6gA4PVyeaPVBV+/H+36epCkOuAxgVzjp25v1Np2lTUga5OTvSJ3f0ohAfN6WLBGCRZgyMIV931ZDra1H7AWBUCD7syD9JGfTe36fF81cndaEuEX5KFwnAYvm6u9CcQaqeL0s2ncYN5wCMCMGHnTifXUzzv08UCaa394uim3pGKF0kAIs3c2AM+bg706krhbTuaLrSxQGwGQg+mqGkpISee+45atu2Lbm7u1N4eDjNnj2bLl682KTl5ebm0vz58yk6Oprc3NzE44IFC+jq1at1vufChQt07733UlRUlHgPl2HmzJmUkvLf8NAl5VV0z9f7RdfBOKcsyt74MfXr10/My+/x8/Oj/v3705IlS6iioqJJZQewBseOHaMpU6ZQSEgIeXh4UJcuXejdd9+l6mr9PVr8PFzUuR/v/nVK3fOlqKiIVqxYQfPmzaO+ffuK/cjBwYFeeOEFs34fAKslQZOUlJRI/fr143pYqWXLltLNN98s9enTR/w/JCRESk5ObtTyMjMzpfj4ePH+uLg4sbxOnTqJ/7dt21bKzs6u9Z4jR45IwcHBYp6YmBjpxhtvlBISEsT/fX19pYMHD4r5nvz5kBT9xBqpx0sbpZffXCRej46OlkaMGCHdcsst4tHd3V1MHzp0qFRWVma09QRgKXbs2CF5eHiI7Zz3Vd7HWrRoIf4/ZcoUqbq6Wu/7rhaXS91e3CD2of/bmSqmJSYmivfp/j3//PNm/lYA1gnBRxM988wz4mDTv39/qaCgQD190aJF6pN4Y0ybNk28jwOIiooK9fR58+aJ6TNmzNCanw+UXbp0Ea/Nnj1b6z3vv/++mN6xY0fptwNp4qAZ8+QaadvpTBEU6QuM0tPTpc6dO4v3LVmypJFrA8CylZeXS7GxsWL7Xrx4sXo677u8D/P0ZcuW1fn+5dtTxH7U/aWNIhg5c+aMNGfOHOnjjz+W9u/fL7300ksIPgAaAcFHE3DNgJ+fnzjYHDhwoNbrXbt2Fa/t27fPoOVdunRJcnR0lFxdXUUQoKm0tFTUpDg5OUlXrlxRT9+6dav4jMDAQK3gRzZgwADxetTU58VB8831Jxosx4oVK8R7Jk2aZFC5AazFDz/8ILZtrhnUxcEDv8bBd13KK6ukEYs2i33p5TXHar3+2muvIfgAaATkfDTB9u3bKS8vj1q3bk3du3ev9frkyZPF4+rVqw1a3vr160Wb8+DBgyksTPtW9tyWfP3111NVVRWtXbtWPX3//v3isWfPnuTt7V1rmUOGDhOP2ce3U/cof1pwbdsGy+Hi4iIeXV1dDSo3gLX4448/tPZNTT169KC4uDg6evQopaam6n2/i5Mj/W+86m7Py3ekUmpWkYlLDGDbEHw0waFDh9QHLX3k6YcPHzbZ8jjhjQUEBOh9z7FsVWJcVVYqvX9Ld3HwbCjZddGiReL5+PHjDSo3gD3ts8PahdLQtiFUUSXR/349Ku4IDQBNg+CjCc6fPy8eIyL0d1eVp587d85ky+Ns/bo+Y0dyFm3ef0w8dynJpshAz1rznD59WvSKmT59Oo0ePVr0ltm7d6/oOTNt2jSDyg1gb/vsixM6iQH6tp3JopX7L5igpAD2wVnpAlijwsJC8ejpWfukzry8vMRjQUGByZY3ZMgQ8cgBw/Hjx6ljx47i/zlF5fTgil1UlLRV/L+sWH/18JUrV+irr77Smvbggw/SwoULydERMSnYFmPtszHBXvTQyLb0+rokevmPEzSkbQiF+bqboMQAtg1nGSvVrl07mjRpksgVmTBhAm3atIny8/Np9qKf6Oiyp0kqUR1E6wokBg0aJKqNKysr6ezZs6LJZfny5dSrV686270BgOjOQbHUpZUf5ZVUiIH7cN8XgMZD8NEEcoJncXGx3tflfAwfHx+TLu+LL74QNSDJyck0YsQIMVjYb8/fThXpZ+ixZxfWmxMic3JyotjYWHr44Ydp2bJlojmGB04CsCXG3GednRzpvVu6kZerE+06m6O+ZQEAGA7BRxNwfoQ8uqg+8nQeodSUy+PAYvPmzSKTf9Z9D5Jvt7HkP2wmLfpxE40c1FvM06lTJ4O/F9ek8EGae9+Ul5cb/D4Ae9tn40K86dUbu6jvepuUnm+0sgLYAwQfTZCQkCAeDxw4oPd1eXrXrl1Nvjwe0nnIiFGUGjeJAkbPpZtm3k/zJ/SlHTt2iNeHDVN1uTUELyswMFA0xXDvFwBbYex9lk3s1opmDogRz/84jLveAjQGgo8mGDhwoGji4OaOgwcP1np95cqV4pHH5zDEmDFjRG7G1q1bKSMjQ+u1srIyMV4IN4+MGzeu1ns5b+N/q47Q2awiaunnTm9NThD3nOEmGR6vY8aMGQZ/L879SEtLI19fXwoODjb4fQCWTu4+Lu+bmhITE8W237lzZ4qJUQUThnr2uo40on0oVdbkfWTklxqpxAC2DcFHE/BJ/YEHHhDP586dq24vZosXLxZjBQwdOlQMAKbpgw8+oPbt29NTTz2lNb1ly5Z06623iqaO+++/X9Q8yB5//HHKzMyk22+/nUJDQ7Xed+rUKfpqywn69eAlcnJ0oPdv7U5SWSFNnTpVdC18+umna3Ut5JvHpafXvjvnyZMn6bbbbhPBDHe/5WAHwFZwkyLnNvF4H++88456Ou+7vA+zRx55pNb7OJeK99k9e/boXa6834X7q3q8/JJ4kY5dyjPZ9wCwGY0ZDhW0byzXt29frRvLyf+v68ZyPPSyvvu0yDeWa926tXidH6dOnaq+10qbNm303ljugUeelBycXSW3iE5Sz+HjpbFjx0re3t7iPTNnzpSqqqpqvYdvKMdDuXfv3l3cTGvy5MlS7969xTR+35AhQ/QO1w5g7bZv366+sRzvq7zP8r7L/+f9QN+N5Xh/4df/+eefWq/dcMMNYjn8FxERKeZz8g6SPFq1lzom9BSvA4B+CD6aobi4WHr22WdFsMD3ZeE7ZPJJPy0tTe/89QUfjAMMvpFcZGSkWB4/Pvjgg1Jubm6teQtLK6Rec9+TPNr0kzwDQsX8fJ+XUaNGSb/88kudZf7666+l2267Tdwpl+986+LiIoWFhUmjR4+Wli9frjdgAbAVR48elW666SYpKChI3MmZ7xzNN5qra7uvL/iQX6vrL6RlhBm+EYB1cuB/lK59gcbhn2zutwdo7ZF0CvFxo7UPDhaPAKCsssoqenzlYfrt4CXx/9kDY+mpce0bvL0BgL3BHmGFPtqcLAIPFycHWjqtBwIPAAvh5uxE79zcjR4YHi/+/+X2FLrts11IRAXQgZoPK7Mp6QrN+Wof8a/26qQudFtf1fgFAGBZNh5Lp0d+PEQFZZXiAuHD23pQn9hApYsFYBEQfFiR/edy6fbPd1NJRRVN6xtFr0xSDXIEAJYpJauI7l2xn05eKRA9Y54c057uHBwrxtQBsGcIPqzEyfQCuvmTneJ+Enxb78+m9yJXZ7SaAVi64vJKevqXI6JLPBvZMYzenpxAfp4uShcNQDEIPqzk6mnqJzspo6CMekYH0Io5fcjTFTckBrAWfJj9evd5Wrj6OJVXVVNEgAd9NK0HdY3wV7poAIpA8GHhjl7Moxlf7qHsonJqF+ZDP97TH1dMAFbqyIU8uv/b/ZSWU0KuTo70v+s60B39otEMA3YHwYcF23Emi+5esZ8KyyqpcytfWj6rDwV7o2cLgDXjptPHfjpEG49fEf8f37UlvX5jF/Jxx0UF2A8EHxaIf5Ivt6fSa2tPiHtGDGgdRJ/c0RMHJwAb2se/2JZCr69LEvt4bLCX6A3TMdxX6aIBmAWCDwtzJb+U/vfrUfqz5qpoQkI4vTm5K7m74F4rALbYg+2Bbw/Q5bxScnN2pJcmdqKbe0WiGQZsHoIPC1FdLdEP+9Lo1T9OiHEBnB0d6H/jO9CMATE4EAHYsJyicnr4x4O0+WSm+P+NPVrRyzd0RlI52DQEHwrj1c+1HIv/PEVJ6QViWkKEH70xuSu1b4EqWAB7ufhYuiWZFm08SdUSUZtQb1p6ew+KD/VRumgAJoHgQyEl5VW0+vAl+mpHKh27lC+m+bg50/xr29CsgbFiQCIAsC+7zmbTvO8SKbOgTDTDzB0eT/cMjRPDtgPYEgQfZlRVLdHe1Bxac/gSrT50WWS9M09XJ5o1MIbuGhxH/p6uShcTABTEgQc3w2w9nSX+HxPkSU+P6yAGJ0MTLNgKBB8mVlRWKa5muD13/bF0cWCR8UBDt/eLFglmgV7WF3RkZqraqAGUFhISQraED8urD1+ml9ccF4MLsq4RfrTg2jY0vF0oghCwegg+TNB2e/xyPv17OpP+PZUpstkrqv5bxb7uzjS6Uwu6LiGcBsUHW3XzCg6AYCls9TBWUFpBSzcn0/IdqVRcXiWmxYV40e19o+mmnhHk54Hu92CdEHw0U2VVtcjZ2JOSQ7tTssVjfmml1jxcwzGkbQhd2yGUBsWH2Mw9WRB8gKWw9cNYdmEZffLvWfp293kx6CDjEVKHtguh67q2pGHtQhGIgFVB8GEgXk1ZheWUlltMJy7ni4Dj2MU80UOlrLJaa17O4egfFyQCDv7jNltbPFHb4ncC62QvhzEOPFYlXqRvdp1T945jXIHaJcJfHHe4eaZLKz9x0YN9FCwVgg+NZNBLV0voQm4JXbxaQhdzS8T/L2r8lesEGZpNKX1iA8Vf39gg6hTuS85OtlG7UR8c2MBS2NthjL/vySsFtObQZVp39DIlZxbVmodrQlqHeFFEgCdFBnpQuL8HBXm5UZC3KwV4ulKQl6uYx9GKm37BetlV8MHBA9dcpGYVUWp2MZ3PLqJzOcV0LruYLuQWa+Vm6MPn2jAfd2oT5k2dW/mJIKNzuB9FBXra5Q6MhFOwFLaWcNpYl/NKaPuZbNqXmiNqZU+mF4i75zaEc8442T3Ux41CfNzUj638PSkm2JPigr0pzNcNFxpgW8FHaUUVVVRVk4uTo9gJeFTPpmzkXGtRXF4pErJyi8tFj5KM/DLKLCyjy1dLKCVbFXBw7QXPWxduQ20V4EGt/FV/fKUg/5+rMMN83W0mXwMAbPtC63RGgfrC6kJNTS6Ppsp/fJfsAp3ctLpwM3JMkBfFh3qLwc/EY5g3RQd5iWM3gNUFH+/9dZre+euU1jSuQHB2/C8YcXJyoLrCEY4jSiqq6mwOqWtH4p0mNtiTogK9KDrIs+bPi1r4ult17xMAAEPxcVO+WBMXbAWlNY9llJZTTClZRZSWW/cFGx+fY4K9KD7kv4CkdYg3RQV5igETUVsC9VH05gFV1bWDBt7ORXWhqldZo3DcwG2YoT7uFOrrRiHebhTi60axQV5iJ+E7R3K1InYKALB3XIvLtbn811BTdXJGIZ3JLKQzGYWq5xmFVFReJR75j47VXjYff4O9XUWzjoerE7k7O5F7zSMfgqslifjSl69/+bgv/i/yWeSlqJ5oXh5z8zbXUPPyXZz4uRO5OGtOk/8cxAixLjrTtac5iEd+rzytqbXvmvj7cMDG5zFef/Ijd0woq+DHKtVz8X+N55VVIlG4a4Q/2QNFaz74B+JmF37k20pXajzXnFYX3kb4bq98Ayau0eANC4EFAIBp8WmD78R7uib4UP0ViMfcYtXIzdbqv4DEgZwcNZuVtE+VmmdOflqhEWw09az60LVtxS027IFdJZwCAIDpc/m4+SarUNWcw007pRXVYrp4rFRVa/NloqMD1zSoes7J/xev1VxDOmg+d3CgyirVBSv/cW2B/JxP+txhQF3boPkaT6/k/EJ+lKc1P1AwFJefgxm+OHZzUV0ku4k/J3Jz0Xju7CgGn5yQEE72AMEHAADYLa5d1wxc/gtmqkXtu2Zlum4GouZr3GQjak1qak7kR84jRI18bQg+AAAAwKzQTwoAAADMCsEHAAAAmBWCDwAAADArBB8AAABgVgg+AAAAwKwQfAAAAIBZIfgAAAAAs0LwAQAAAGaF4AMAAADMCsEHAAAAmBWCDwAAADArBB8AAABgVgg+AAAAwKycDZmJb3xbXl5u+tIAAACA0bm6upKDgwNZVfDBgcfrr79u+tIAAACA0c2dO5eCg4PJqoIPjpiefPJJMrf09HRavnw5zZw5k1q0aEH2DOtCBevhP1gXKlgP/8G6UMF6qL0uLK31wqDgg6tq3NzcyNw46JEflfh8S4J1oYL18B+sCxWsh/9gXahgPdReF5bU5MKQcAoAAABmZdHBh7e3Nw0dOlQ82jusCxWsh/9gXahgPfwH60IF68Hy14WDxF1ZAAAAAMzEoms+AAAAwPYg+AAAAACzQvABAAAAZoXgAwAAAMwKwQcAAABYV/Cxd+9eGjduHPn7+5OXlxf169ePfvzxx0Yto6ysjF566SVq06YNubu7U3h4ON19992UkZGhd/6SkhJavHgx9ejRgwICAsRnJyQk0CuvvEJ5eXm15h82bJgYYEXfX0xMDBmDEushNzeXHn30UYqPjxcD6YSEhNDkyZPp2LFjdX7GqVOn6OabbxbD7Hp4eIj1tnTpUnH/HmOxhnVhDdtEcnIyvfDCCzRhwgRq1aqVwWXbsGGD6Frn4+NDvr6+NHz4cPr777+tepswx7qw1W1iyZIlNGvWLOratSs5OzuL92zevLne91y+fJnmzJlDLVu2FPtfu3btxPG1oqKCjMUa1gWPkFrXNmGsQbv2mnk9nD59ml599VUaMmSIOK7yIGSRkZE0ffp0SkpKMt82ITXDpk2bJBcXF8nHx0e66667pIcffliKjo7mI5b09ttvG7SMqqoqafTo0eI9/fr1k5544gnpxhtvlBwcHKS4uDgpIyNDa/7y8nKpb9++Yv5u3bpJCxYsEH8JCQliWqdOnaSioiKt9wwdOlS89vzzz9f6e+edd5qzChRbD1lZWVKbNm3E/P379xefeeutt0qurq6Sp6entGvXrlqfcezYMcnPz0/Mc/vtt0uPP/64WF+8jAceeKDZ68Ga1oU1bBPLli0T8zs5OUmdO3eWHB0dxTLqs2LFCvGekJAQ8ZvyHz/ndffTTz9Z7TZhjnVhq9sEz89/LVu2lFq0aCGe//PPP3XOf/nyZSkyMlKsJ97veP/j/ZDfN2HCBKm6ulqyl3UxY8YMMc/8+fP1bhfWuB6mTp0q5ud57733XrHPjx07Vkzz8PCQtmzZYpZtosnBR0VFhdS6dWvJzc1NSkxMVE+/evWq1LZtW3EwS01NbXA5X375pfgCfLLQ/AJLly4V0++++26t+X/44QcxfdKkSbWWNXHiRPHaV199pfegYgpKrYe5c+eK6byxatqxY4fYCDt27ChO4pqGDBki3rN27Vr1tLKyMmnw4MFiOr/XXtaFNWwTycnJ0s6dO6Xi4mLxf15efQeVnJwcyd/fXwoODpbS0tLU0/k5T+O//Px8q9wmzLEubHGbYGvWrBEnD3bPPfc0eMKdPn26mIf3Nxnvh7fccouY/u2330r2si7k4CMlJUUytgqF1gMHKwcOHKg1/bvvvhPflY+X5tgmmrynbdiwQXzorFmzar22fPly8dqLL77Y4HL4SpXn1V3J/MX4KtfLy0u9Utlrr70m5v/0009rLYun6YsYTXlQUWo9REREiAi3oKCg1rJuuOEGsSyOqmUnT54U04YPH15r/s2bN9f5HWxxXVjLNqGroYPKJ598UueyX3jhhVqBuTVtE6ZeF7a6Tehq6ITLARkvk/cz3atZ3h/r2l5scV2YOvjYYCHrQRMHPfy5mZmZJt8mmpzzIbeTjRo1qtZro0ePFo9btmypdxmlpaW0e/du0XYUHR2t9Rq3W40cOZKKiopo37596umdO3cWj+vWrau1vD/++EO8j9t19fn2229FW9e7774ryl9dXU3NpdR64DsVchu9viFzY2NjxeOmTZsMKuegQYNEW2ND5bSVdWEN24Q5Ptdatglzf64tbRONtXPnTpFvxfuZbk4D74+8X27fvp2qqqpsfl1oWrNmDb322msi15DPPca4Q+xmC1wPLi4u4pHzYUy9TRh0V9u6klYYJwTq4lsY84lAnqe+RBnesfUtQ3PZvJzBgweL5+PHj6cbbriBVq1aRd27dxdJYuyff/6hlJQU+vTTT0Uiqj7Tpk3T+n/btm3pm2++oV69epG1rQc+2XLyZWFhYa2TLq8HOZHQkHI6OTmJk/Tx48epsrJSa8OzxXVhDduEsT9Xc71Z2zZhjnVhq9uEMcspTz958iSdO3eO4uLibHpdaJo3b57W/znpctmyZeogwRbWw549e0SCfu/evUXyqyHlbM420eSaD7lXiZ+fn97XObNcX8+Txi5Dcz7GkdfPP/9MTzzxBB06dEhcnfAfP580aZKIznRNnDhRRK4XL16k4uJicUCdP3++ONHx/OfPn2/EN7eM9TB27Fhxkn7xxRe15uVaA/6u7OrVq436DF5eQUFBvWW1hXVhDduEsT9X33qzlm3CHOvCVreJxmrK/mer64JxjxDuecK/Pfey5BMx98Lj4wn3LtGsgbXm9ZCXl0czZswgR0dHevPNNxtdTs35bHacDz4ocJCxfPly+u677ygrK0v8ff/997R+/Xrq06cPpaamar3noYceEjUm3K2IuxJ26NBBBCxPP/202Ijefvttsja8A3D0zWXnKnLuZspXbLyzdOzYUczDG5I9aMq6sMVtApoH2wTomj17Nk2ZMkV0ReXupdyV/9lnn6UPPvhANL3wscfalZSUiHMqd7NduHChujXB1Jp8dpKjoLqinfz8/DojpcYsQ3M+xm2xv//+u2hemTp1KgUFBYk/fv7JJ5+I6nfue2yIe+65Rzxye5W1rYeIiAjRP5z7XXPTwvvvv0+7du0SOwMfLFloaGijPoNrlXg8BFtfF9awTRj7c/WtN2vZJsyxLmx1m2ispux/trou6sO1BByMWPs2UVpaKmr8OG3hqaeeUh8vG1tOzflMHnzU127KCYDc/l5XG5GM24f4irSudi19bU1yoqm+pFJ5WmJiokHfgYMWPrhyAqO1rQfGA8p8/vnnopqYo3CuHubmqBMnTojXNduo6ysnJwrxSZvb+Jvatm9N68Iatgljf66+9WYt24SxP7ehNmxb2iYaq6F8GJ7Og1JFRUXZ/LqoD+dEcV6ENW8TJSUlounozz//pMcff1xc2De2nM3ZJpocfPCogWzjxo16RxXUnKcuXLXJzSRysoom7gbMK4Uz7jVPHHKWcWZmZq3lydN4hEtDcIINf05zRi9Uaj3UhU8a3ATFJ4ybbrrJoHJu27ZN7EQNldNW1oU1bBPm+Fxr2SaU/lxr3iYai0fX5BMJ72e6I9zy/sj75cCBA5sVkFrLuqgP54BwgGCt20RJSYmo8eDfmZup33jjDfNvE80ZIIX7/dY3QIpm3+hLly5JJ06cEK83Z0ApuW82D3qiOXBUZWWlNG3aNPHaM888o55+9uxZKTs7u1b5L1y4oB7JUbe/vzWsBx7pVXOsC8brg0d75fkfeuihWmVtaECp7du3N3k9WNO6sJZtQpchA2vxaKXGHGTMUrYJU68LW90mdDVnkDHeH401yJg1rAsejIx/f125ubliXIumjsOh9HooKSmRRo4cqXdgRnNuE2YbXl0erIVHV2toKO2bbrpJDOMaGxtbayjtc+fOqYfF5YPCvHnzxB+PysbTeJhtPvDI+PPc3d3FyuaTFi//5ptvFgNV8fwcsDR3uGAl1gMfRH19faXJkydLjz32mBj+t3379uL948ePl0pLS2uV8+jRo+qhtO+44w7Fh9JWal1YyzbBA/3wa/IfD6TGZdScpjkYUENDiv/4449Wu02Yel3Y8jbBAzPKr8mDSPH+JU9btWqV1vx8kpOH0ub9TnMo7euvv97sw6srtS44KHF2dhYB+pw5c6Qnn3xSbAdBQUHifddcc404kVvbephRsxw+j+obMp7/dAdVM8U20ezh/Hbv3i2NGTNGHPx5XPg+ffpI33//fa356lpxjE8OPOogDzXLB0FeKXfeeaeUnp6u9zMvXrwoDiTx8fFifo702rVrJ048moEHO3TokDiocnDCwy3zxsRXPaNGjdJbTmtZD3zVxt+LI2c+aPLGyyODfvbZZ7WGEteUlJQkTtKBgYFivXXp0kX68MMPjXJAsZZ1YS3bBB8AeHp9f/pGXly3bp2oueADkLe3txi5888//7TqbcLU68KWtwl55Na6/vTdo4RPNrNnz5bCwsLE/scXdQsXLhS1YvayLs6fPy+OOXzfMA44eJvgbYODkY8//ljUtlvjehjawDqoqzbI2NuEA//T1HYjAAAAgMayj4EgAAAAwGIg+AAAAACzQvABAAAAZoXgAwAAAMwKwQcAAACYFYIPAAAAMCsEHwAAAGBWCD4AAADArBB8AAAAgFkh+AAAAACzQvABAAAAZoXgAwAAAMic/h+ta1+spwSQIAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "az.plot_posterior(trace, var_names=[\"delta\"], hdi_prob=0.90)\n",
    "plt.title(f\"Posterior of {metric} difference: {model_2} - {model_1}\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "3c2c5a42",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(1.0)"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "delta_samples = trace.posterior[\"delta\"].values.flatten()\n",
    "(delta_samples > 0.00).mean()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pymc_env",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
